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

Diff of /MITgcm/model/src/external_forcing_surf.F

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

revision 1.16 by edhill, Fri Oct 24 05:29:35 2003 UTC revision 1.49 by jmc, Mon Dec 21 00:24:58 2009 UTC
# Line 3  C $Name$ Line 3  C $Name$
3    
4  #include "PACKAGES_CONFIG.h"  #include "PACKAGES_CONFIG.h"
5  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
6    
7  CBOP  CBOP
8  C     !ROUTINE: EXTERNAL_FORCING_SURF  C     !ROUTINE: EXTERNAL_FORCING_SURF
9  C     !INTERFACE:  C     !INTERFACE:
10        SUBROUTINE EXTERNAL_FORCING_SURF(        SUBROUTINE EXTERNAL_FORCING_SURF(
11       I             bi, bj, iMin, iMax, jMin, jMax,       I             bi, bj, iMin, iMax, jMin, jMax,
12       I             myTime, myIter, myThid )       I             myTime, myIter, myThid )
13  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
14  C     *==========================================================*  C     *==========================================================*
15  C     | SUBROUTINE EXTERNAL_FORCING_SURF                            C     | SUBROUTINE EXTERNAL_FORCING_SURF
16  C     | o Determines forcing terms based on external fields        C     | o Determines forcing terms based on external fields
17  C     |   relaxation terms etc.                                    C     |   relaxation terms etc.
18  C     *==========================================================*  C     *==========================================================*
19  C     \ev  C     \ev
20    
# Line 28  C     === Global variables === Line 28  C     === Global variables ===
28  #include "DYNVARS.h"  #include "DYNVARS.h"
29  #include "GRID.h"  #include "GRID.h"
30  #include "SURFACE.h"  #include "SURFACE.h"
31    #ifdef ALLOW_SEAICE
32    #include "SEAICE_PARAMS.h"
33    #include "SEAICE.h"
34    #endif /* ALLOW_SEAICE */
35    #ifdef ALLOW_SHELFICE
36    #include "SHELFICE.h"
37    #endif /* ALLOW_SHELFICE */
38    
39  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
40  C     === Routine arguments ===  C     === Routine arguments ===
41  C     myTime - Current time in simulation  C     bi,bj  :: tile indices
42  C     myIter - Current iteration number in simulation  C     iMin,iMax, jMin,jMax :: Range of points for calculation
43    C     myTime :: Current time in simulation
44    C     myIter :: Current iteration number in simulation
45  C     myThid :: Thread no. that called this routine.  C     myThid :: Thread no. that called this routine.
       _RL myTime  
       INTEGER myIter  
       INTEGER myThid  
46        INTEGER bi,bj        INTEGER bi,bj
47        INTEGER iMin, iMax        INTEGER iMin, iMax
48        INTEGER jMin, jMax        INTEGER jMin, jMax
49          _RL myTime
50          INTEGER myIter
51          INTEGER myThid
52    
53  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
54  C     === Local variables ===  C     === Local variables ===
55    C     i,j    :: loop indices
56    C     ks     :: index of surface interface layer
57        INTEGER i,j        INTEGER i,j
58  C     number of surface interface layer        INTEGER ks
       INTEGER kSurface  
59  CEOP  CEOP
60    #ifdef ALLOW_DIAGNOSTICS
61          _RL tmpFac
62    #endif /* ALLOW_DIAGNOSTICS */
63    
64        if ( buoyancyRelation .eq. 'OCEANICP' ) then        IF ( usingPCoords ) THEN
65         kSurface        = Nr         ks        = Nr
66        else        ELSE
67         kSurface        = 1         ks        = 1
68        endif        ENDIF
69    
70  C--   Surface Fluxes :  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
71    
72        DO j = jMin, jMax        IF ( doThetaClimRelax .OR. doSaltClimRelax ) THEN
73           DO i = iMin, iMax  C--   Start with surface restoring term :
74    
75  c     Zonal wind stress fu:  #ifdef ALLOW_SEAICE
76            surfaceTendencyU(i,j,bi,bj) =         IF ( useSEAICE .AND. (.NOT. SEAICErestoreUnderIce) ) THEN
77       &      fu(i,j,bi,bj)*horiVertRatio*recip_rhoConst  C     Do not restore under sea-ice
78       &           *recip_drF(kSurface)*recip_hFacW(i,j,kSurface,bi,bj)          DO j = jMin, jMax
79  c     Meridional wind stress fv:           DO i = iMin, iMax
80            surfaceTendencyV(i,j,bi,bj) =  C     Heat Flux (restoring term) :
81       &      fv(i,j,bi,bj)*horiVertRatio*recip_rhoConst            surfaceForcingT(i,j,bi,bj) =
82       &           *recip_drF(kSurface)*recip_hFacS(i,j,kSurface,bi,bj)       &      -lambdaThetaClimRelax(i,j,bi,bj)*(1.-AREA(i,j,bi,bj))
83  c     Net heat flux Qnet:       &         *(theta(i,j,ks,bi,bj)-SST(i,j,bi,bj))
84            surfaceTendencyT(i,j,bi,bj) =       &         *drF(ks)*_hFacC(i,j,ks,bi,bj)
85       &      -Qnet(i,j,bi,bj)*recip_Cp*horiVertRatio*recip_rhoConst  C     Salt Flux (restoring term) :
86       &           *recip_drF(kSurface)*recip_hFacC(i,j,kSurface,bi,bj)            surfaceForcingS(i,j,bi,bj) =
87  C     Net Salt Flux :       &      -lambdaSaltClimRelax(i,j,bi,bj) *(1.-AREA(i,j,bi,bj))
88            surfaceTendencyS(i,j,bi,bj) =       &         *(salt(i,j,ks,bi,bj)-SSS(i,j,bi,bj))
89       &      -saltFlux(i,j,bi,bj)*horiVertRatio*recip_rhoConst       &         *drF(ks)*_hFacC(i,j,ks,bi,bj)
90       &           *recip_drF(kSurface)*recip_hFacC(i,j,kSurface,bi,bj)           ENDDO
91            ENDDO
92  #ifdef ALLOW_PASSIVE_TRACER         ELSE
93  c ***  define the tracer surface tendency here ***  #endif /* ALLOW_SEAICE */
94  #endif /* ALLOW_PASSIVE_TRACER */          DO j = jMin, jMax
95             DO i = iMin, iMax
96    C     Heat Flux (restoring term) :
97              surfaceForcingT(i,j,bi,bj) =
98         &      -lambdaThetaClimRelax(i,j,bi,bj)
99         &         *(theta(i,j,ks,bi,bj)-SST(i,j,bi,bj))
100         &         *drF(ks)*_hFacC(i,j,ks,bi,bj)
101    C     Salt Flux (restoring term) :
102              surfaceForcingS(i,j,bi,bj) =
103         &      -lambdaSaltClimRelax(i,j,bi,bj)
104         &         *(salt(i,j,ks,bi,bj)-SSS(i,j,bi,bj))
105         &         *drF(ks)*_hFacC(i,j,ks,bi,bj)
106             ENDDO
107            ENDDO
108    #ifdef ALLOW_SEAICE
109           ENDIF
110    #endif /* ALLOW_SEAICE */
111    
112    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
113    #ifdef NONLIN_FRSURF
114    C-    T,S surface forcing will be applied (thermodynamics) after the update
115    C     of surf.thickness (hFac): account for change in surf.thickness
116           IF (staggerTimeStep.AND.nonlinFreeSurf.GT.0) THEN
117            IF (select_rStar.GT.0) THEN
118    # ifndef DISABLE_RSTAR_CODE
119             DO j=jMin,jMax
120              DO i=iMin,iMax
121                surfaceForcingT(i,j,bi,bj) = surfaceForcingT(i,j,bi,bj)
122         &                                  * rStarExpC(i,j,bi,bj)
123                surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
124         &                                  * rStarExpC(i,j,bi,bj)
125              ENDDO
126           ENDDO           ENDDO
127        ENDDO  # endif /* DISABLE_RSTAR_CODE */
128                    ELSE
129  C--   Surface restoring term :           DO j=jMin,jMax
130              DO i=iMin,iMax
131               IF (ks.EQ.ksurfC(i,j,bi,bj)) THEN
132                surfaceForcingT(i,j,bi,bj) = surfaceForcingT(i,j,bi,bj)
133         &             *_recip_hFacC(i,j,ks,bi,bj)*hFac_surfC(i,j,bi,bj)
134                surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
135         &             *_recip_hFacC(i,j,ks,bi,bj)*hFac_surfC(i,j,bi,bj)
136               ENDIF
137              ENDDO
138             ENDDO
139            ENDIF
140           ENDIF
141    #endif /* NONLIN_FRSURF */
142    
143    #ifdef ALLOW_DIAGNOSTICS
144           IF ( useDiagnostics ) THEN
145    
146    C     tRelax (temperature relaxation [W/m2], positive <-> increasing Theta)
147            tmpFac = HeatCapacity_Cp*rUnit2mass
148            CALL DIAGNOSTICS_SCALE_FILL(
149         &           surfaceForcingT(1-OLx,1-OLy,bi,bj),tmpFac,1,
150         &                             'TRELAX  ',0, 1,2,bi,bj,myThid)
151    
152    C     sRelax (salt relaxation [g/m2/s], positive <-> increasing Salt)
153            tmpFac = rUnit2mass
154            CALL DIAGNOSTICS_SCALE_FILL(
155         &           surfaceForcingS(1-OLx,1-OLy,bi,bj),tmpFac,1,
156         &                             'SRELAX  ',0, 1,2,bi,bj,myThid)
157    
158           ENDIF
159    #endif /* ALLOW_DIAGNOSTICS */
160    
161          ELSE
162    C--   No restoring for T & S : set surfaceForcingT,S to zero :
163    
       IF ( doThetaClimRelax .OR. doSaltClimRelax ) THEN  
164         DO j = jMin, jMax         DO j = jMin, jMax
165          DO i = iMin, iMax          DO i = iMin, iMax
166  C     Heat Flux (restoring term) :            surfaceForcingT(i,j,bi,bj) = 0. _d 0
167           IF ( abs(yC(i,j,bi,bj)).LE.latBandClimRelax ) THEN            surfaceForcingS(i,j,bi,bj) = 0. _d 0
           surfaceTendencyT(i,j,bi,bj) = surfaceTendencyT(i,j,bi,bj)  
      &      -lambdaThetaClimRelax  
      &         *(theta(i,j,kSurface,bi,bj)-SST(i,j,bi,bj))  
 C     Salt Flux (restoring term) :  
           surfaceTendencyS(i,j,bi,bj) = surfaceTendencyS(i,j,bi,bj)  
      &      -lambdaSaltClimRelax  
      &         *(salt(i,j,kSurface,bi,bj)-SSS(i,j,bi,bj))  
          ENDIF  
168          ENDDO          ENDDO
169         ENDDO         ENDDO
170    
171    C--   end restoring / no restoring block.
172        ENDIF        ENDIF
173    
174  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
 C--   Fresh-water flux  
175    
176  #ifdef NONLIN_FRSURF  C--   Surface Fluxes :
177        IF ( (nonlinFreeSurf.GT.0 .OR. buoyancyRelation.EQ.'OCEANICP')  
178          DO j = jMin, jMax
179             DO i = iMin, iMax
180    
181    C     Zonal wind stress fu:
182              surfaceForcingU(i,j,bi,bj) =
183         &      fu(i,j,bi,bj)*mass2rUnit
184    C     Meridional wind stress fv:
185              surfaceForcingV(i,j,bi,bj) =
186         &      fv(i,j,bi,bj)*mass2rUnit
187    C     Net heat flux Qnet:
188              surfaceForcingT(i,j,bi,bj) = surfaceForcingT(i,j,bi,bj)
189         &       - ( Qnet(i,j,bi,bj)
190    #ifdef SHORTWAVE_HEATING
191         &          -Qsw(i,j,bi,bj)
192    #endif
193         &         ) *recip_Cp*mass2rUnit
194    C     Net Salt Flux :
195              surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
196         &      -saltFlux(i,j,bi,bj)*mass2rUnit
197    
198             ENDDO
199          ENDDO
200    
201    #ifdef ALLOW_SALT_PLUME
202    C saltPlume is the amount of salt rejected by ice while freezing;
203    C it is here subtracted from surfaceForcingS and will be redistributed
204    C to multiple vertical levels later on as per Duffy et al. (GRL 1999)
205          IF ( useSALT_PLUME ) THEN
206             CALL SALT_PLUME_FORCING_SURF(
207         I        bi, bj, iMin, iMax, jMin, jMax,
208         I        myTime,myIter,myThid )
209          ENDIF
210    #endif /* ALLOW_SALT_PLUME */
211    
212    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
213    C--   Fresh-water flux
214    
215    C-    Apply mask on Fresh-Water flux
216    C      (needed for SSH forcing, whether or not exactConserv is used)
217          IF ( useRealFreshWaterFlux ) THEN
218            DO j=1-OLy,sNy+OLy
219             DO i=1-OLx,sNx+OLx
220               EmPmR(i,j,bi,bj) = EmPmR(i,j,bi,bj)*maskInC(i,j,bi,bj)
221             ENDDO
222            ENDDO
223          ENDIF
224    
225    #ifdef EXACT_CONSERV
226    C NB: synchronous time step: PmEpR lag 1 time step behind EmPmR
227    C     to stay consitent with volume change (=d/dt etaH).
228          IF ( staggerTimeStep ) THEN
229            DO j=1-OLy,sNy+OLy
230             DO i=1-OLx,sNx+OLx
231               PmEpR(i,j,bi,bj) = -EmPmR(i,j,bi,bj)
232             ENDDO
233            ENDDO
234          ENDIF
235    
236          IF ( (nonlinFreeSurf.GT.0 .OR. usingPCoords)
237       &     .AND. useRealFreshWaterFlux ) THEN       &     .AND. useRealFreshWaterFlux ) THEN
238    
239  c-  NonLin_FrSurf and RealFreshWaterFlux : PmEpR effectively changes  C--   NonLin_FrSurf and RealFreshWaterFlux : PmEpR effectively changes
240  c   the water column height ; temp., salt, (tracer) flux associated  C     the water column height ; temp., salt, (tracer) flux associated
241  c   with this input/output of water is added here to the surface tendency.  C     with this input/output of water is added here to the surface tendency.
 c  
 c NB: PmEpR lag 1 time step behind EmPmR ( PmEpR_n = - EmPmR_n-1 ) to stay  
 c     consitent with volume change (=d/dt etaN).  
242    
243         IF (temp_EvPrRn.NE.UNSET_RL) THEN         IF (temp_EvPrRn.NE.UNSET_RL) THEN
244          DO j = jMin, jMax          DO j = jMin, jMax
245           DO i = iMin, iMax           DO i = iMin, iMax
246            surfaceTendencyT(i,j,bi,bj) = surfaceTendencyT(i,j,bi,bj)            surfaceForcingT(i,j,bi,bj) = surfaceForcingT(i,j,bi,bj)
247       &      + PmEpR(i,j,bi,bj)       &      + PmEpR(i,j,bi,bj)
248       &         *( temp_EvPrRn - theta(i,j,kSurface,bi,bj) )       &          *( temp_EvPrRn - theta(i,j,ks,bi,bj) )
249       &           *recip_drF(kSurface)*recip_hFacC(i,j,kSurface,bi,bj)       &          *mass2rUnit
      &         *convertEmP2rUnit  
250           ENDDO           ENDDO
251          ENDDO          ENDDO
252         ENDIF         ENDIF
# Line 131  c     consitent with volume change (=d/d Line 254  c     consitent with volume change (=d/d
254         IF (salt_EvPrRn.NE.UNSET_RL) THEN         IF (salt_EvPrRn.NE.UNSET_RL) THEN
255          DO j = jMin, jMax          DO j = jMin, jMax
256           DO i = iMin, iMax           DO i = iMin, iMax
257            surfaceTendencyS(i,j,bi,bj) = surfaceTendencyS(i,j,bi,bj)            surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
258       &      + PmEpR(i,j,bi,bj)       &      + PmEpR(i,j,bi,bj)
259       &         *( salt_EvPrRn - salt(i,j,kSurface,bi,bj) )       &          *( salt_EvPrRn - salt(i,j,ks,bi,bj) )
260       &           *recip_drF(kSurface)*recip_hFacC(i,j,kSurface,bi,bj)       &          *mass2rUnit
      &         *convertEmP2rUnit  
261           ENDDO           ENDDO
262          ENDDO          ENDDO
263         ENDIF         ENDIF
264    
 #ifdef ALLOW_PASSIVE_TRACER  
 c  *** add the tracer flux associated with P-E+R here ***  
 c      IF (trac_EvPrRn.NE.UNSET_RL) THEN  
 c    &      + PmEpR(i,j,bi,bj)*( trac_EvPrRn - tr1(i,j,kSurface,bi,bj) )  
 c    &           *recip_drF(kSurface)*recip_hFacC(i,j,kSurface,bi,bj)  
 c      ENDIF  
 #endif /* ALLOW_PASSIVE_TRACER */  
   
265  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
266        ELSE        ELSE
267  #else /* NONLIN_FRSURF */  #else /* EXACT_CONSERV */
268        IF (.TRUE.) THEN        IF (.TRUE.) THEN
269  #endif /* NONLIN_FRSURF */  #endif /* EXACT_CONSERV */
270    
271  c- EmPmR does not really affect the water column height (for tracer budget)  C--   EmPmR does not really affect the water column height (for tracer budget)
272  c   and is converted to a salt tendency.  C     and is converted to a salt tendency.
273    
274         IF (convertFW2Salt .EQ. -1.) THEN         IF (convertFW2Salt .EQ. -1.) THEN
275  c- converts EmPmR to salinity tendency using surface local salinity  C-    use local surface tracer field to calculate forcing term:
276          DO j = jMin, jMax  
277           DO i = iMin, iMax          IF (temp_EvPrRn.NE.UNSET_RL) THEN
278            surfaceTendencyS(i,j,bi,bj) = surfaceTendencyS(i,j,bi,bj)  C     account for Rain/Evap heat content (temp_EvPrRn) using local SST
279       &      + EmPmR(i,j,bi,bj)*salt(i,j,kSurface,bi,bj)           DO j = jMin, jMax
280       &           *recip_drF(kSurface)*recip_hFacC(i,j,kSurface,bi,bj)            DO i = iMin, iMax
281       &         *convertEmP2rUnit             surfaceForcingT(i,j,bi,bj) = surfaceForcingT(i,j,bi,bj)
282         &       + EmPmR(i,j,bi,bj)
283         &           *( theta(i,j,ks,bi,bj) - temp_EvPrRn )
284         &           *mass2rUnit
285              ENDDO
286           ENDDO           ENDDO
287          ENDDO          ENDIF
288         ELSE          IF (salt_EvPrRn.NE.UNSET_RL) THEN
289  c- converts EmPmR to virtual salt flux using uniform salinity (default=35)  C     converts EmPmR to salinity tendency using surface local salinity
290          DO j = jMin, jMax           DO j = jMin, jMax
291           DO i = iMin, iMax            DO i = iMin, iMax
292            surfaceTendencyS(i,j,bi,bj) = surfaceTendencyS(i,j,bi,bj)             surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
293       &      + EmPmR(i,j,bi,bj)*convertFW2Salt       &       + EmPmR(i,j,bi,bj)
294       &           *recip_drF(kSurface)*recip_hFacC(i,j,kSurface,bi,bj)       &           *( salt(i,j,ks,bi,bj) - salt_EvPrRn )
295       &         *convertEmP2rUnit       &           *mass2rUnit
296              ENDDO
297           ENDDO           ENDDO
298          ENDDO          ENDIF
299    
300           ELSE
301    C-    use uniform tracer value to calculate forcing term:
302    
303            IF (temp_EvPrRn.NE.UNSET_RL) THEN
304    C     account for Rain/Evap heat content (temp_EvPrRn) assuming uniform SST (=tRef)
305             DO j = jMin, jMax
306              DO i = iMin, iMax
307               surfaceForcingT(i,j,bi,bj) = surfaceForcingT(i,j,bi,bj)
308         &       + EmPmR(i,j,bi,bj)
309         &           *( tRef(ks) - temp_EvPrRn )
310         &           *mass2rUnit
311              ENDDO
312             ENDDO
313            ENDIF
314            IF (salt_EvPrRn.NE.UNSET_RL) THEN
315    C     converts EmPmR to virtual salt flux using uniform salinity (default=35)
316             DO j = jMin, jMax
317              DO i = iMin, iMax
318               surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
319         &       + EmPmR(i,j,bi,bj)
320         &           *( convertFW2Salt - salt_EvPrRn )
321         &           *mass2rUnit
322              ENDDO
323             ENDDO
324            ENDIF
325    
326    C-    end local-surface-tracer / uniform-value distinction
327         ENDIF         ENDIF
328    
329        ENDIF        ENDIF
330    
331  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
332    
333  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
334        IF ( usePTRACERS ) THEN        IF ( usePTRACERS ) THEN
335           CALL PTRACERS_FORCING_SURF(           CALL PTRACERS_FORCING_SURF(
336       I                              bi, bj, iMin, iMax, jMin, jMax,       I        bi, bj, iMin, iMax, jMin, jMax,
337       I                              myTime,myIter,myThid )       I        myTime,myIter,myThid )
338        ENDIF        ENDIF
339  #endif /* ALLOW_PTRACERS */  #endif /* ALLOW_PTRACERS */
340    
341  #ifdef ATMOSPHERIC_LOADING  #ifdef ATMOSPHERIC_LOADING
342    C-- Atmospheric surface Pressure loading : added to phi0surf when using Z-coord;
343    C   Not yet implemented for Ocean in P: would need to be applied to the other end
344    C   of the column, as a vertical velocity (omega); (meaningless for Atmos in P).
345    C- Note:
346    C   Using P-coord., a hack (now directly applied from S/R INI_FORCING)
347    C   is sometime used to read phi0surf from a file (pLoadFile) instead
348    C   of computing it from bathymetry & density ref. profile.
349    
350  C-- Atmospheric surface Pressure loading :        IF ( usingZCoords ) THEN
351    C   The true atmospheric P-loading is not yet implemented for P-coord
352        IF (buoyancyRelation .eq. 'OCEANIC' ) THEN  C   (requires time varying dP(Nr) like dP(k-bottom) with NonLin FS).
353           IF ( useRealFreshWaterFlux ) THEN
354            DO j = jMin, jMax
355             DO i = iMin, iMax
356              phi0surf(i,j,bi,bj) = ( pLoad(i,j,bi,bj)
357         &                        +sIceLoad(i,j,bi,bj)*gravity
358         &                          )*recip_rhoConst
359             ENDDO
360            ENDDO
361           ELSE
362          DO j = jMin, jMax          DO j = jMin, jMax
363           DO i = iMin, iMax           DO i = iMin, iMax
364            phi0surf(i,j,bi,bj) = pload(i,j,bi,bj)*recip_rhoConst            phi0surf(i,j,bi,bj) = pLoad(i,j,bi,bj)*recip_rhoConst
365           ENDDO           ENDDO
366          ENDDO          ENDDO
367        ELSEIF ( buoyancyRelation .eq. 'OCEANICP' ) THEN         ENDIF
368  C-- This is a hack used to read phi0surf from a file (ploadFile)  c     ELSEIF ( usingPCoords ) THEN
369    C-- This is a hack used to read phi0surf from a file (pLoadFile)
370  C   instead of computing it from bathymetry & density ref. profile.  C   instead of computing it from bathymetry & density ref. profile.
371  C   The true atmospheric P-loading is not yet implemented for P-coord  C   ==> now done only once, in S/R INI_FORCING
372  C   (requires time varying dP(Nr) like dP(k-bottom) with NonLin FS).  c       DO j = jMin, jMax
373    c        DO i = iMin, iMax
374    c         phi0surf(i,j,bi,bj) = pLoad(i,j,bi,bj)
375    c        ENDDO
376    c       ENDDO
377          ENDIF
378    #endif /* ATMOSPHERIC_LOADING */
379    
380    #ifdef ALLOW_SHELFICE
381          IF ( usingZCoords ) THEN
382           IF ( useSHELFICE) THEN
383          DO j = jMin, jMax          DO j = jMin, jMax
384           DO i = iMin, iMax           DO i = iMin, iMax
385            phi0surf(i,j,bi,bj) = pload(i,j,bi,bj)            phi0surf(i,j,bi,bj) = phi0surf(i,j,bi,bj)
386         &         + shelficeLoadAnomaly(i,j,bi,bj)*recip_rhoConst
387           ENDDO           ENDDO
388          ENDDO          ENDDO
389           ENDIF
390        ENDIF        ENDIF
391    #endif /* ALLOW_SHELFICE */
392    
393  #endif /* ATMOSPHERIC_LOADING */  #ifdef ALLOW_EBM
394    c--    Values for surfaceForcingT, surfaceForcingS
395    c      are overwritten by those produced by EBM
396    cph    AD recomputation problems if these IF useEBM are used
397    cph      IF ( useEBM ) THEN
398           CALL EBM_FORCING_SURF(
399         I        bi, bj, iMin, iMax, jMin, jMax,
400         I        myTime,myIter,myThid )
401    cph      ENDIF
402    #endif
403    
404        RETURN        RETURN
405        END        END
Legend:
Removed from v.1.16  
changed lines
  Added in v.1.49
  ViewVC Help
Powered by ViewVC 1.1.22