/[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.9 by mlosch, Wed Oct 2 17:59:57 2002 UTC revision 1.36 by jmc, Sat Jul 22 03:53:13 2006 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    
4    #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             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 26  C     === Global variables === Line 27  C     === Global variables ===
27  #include "FFIELDS.h"  #include "FFIELDS.h"
28  #include "DYNVARS.h"  #include "DYNVARS.h"
29  #include "GRID.h"  #include "GRID.h"
 #ifdef NONLIN_FRSURF  
30  #include "SURFACE.h"  #include "SURFACE.h"
31  #endif  #ifdef ALLOW_SEAICE
32    #include "SEAICE.h"
33    #endif /* ALLOW_SEAICE */
34    #ifdef ALLOW_SHELFICE
35    #include "SHELFICE.h"
36    #endif /* ALLOW_SHELFICE */
37  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
38  C     === Routine arguments ===  C     === Routine arguments ===
39    C     bi,bj  :: tile indices
40    C     iMin,iMax, jMin,jMax :: Range of points for calculation
41    C     myTime :: Current time in simulation
42    C     myIter :: Current iteration number in simulation
43  C     myThid :: Thread no. that called this routine.  C     myThid :: Thread no. that called this routine.
       INTEGER myThid  
44        INTEGER bi,bj        INTEGER bi,bj
45        INTEGER iMin, iMax        INTEGER iMin, iMax
46        INTEGER jMin, jMax        INTEGER jMin, jMax
47          _RL myTime
48          INTEGER myIter
49          INTEGER myThid
50    
51  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
52  C     === Local variables ===  C     === Local variables ===
53    C     i,j    :: loop indices
54    C     ks     :: index of surface interface layer
55        INTEGER i,j        INTEGER i,j
56  C     number of surface interface layer        INTEGER ks
       INTEGER kSurface  
       _RL convertVol2Mass  
57  CEOP  CEOP
58    #ifdef ALLOW_DIAGNOSTICS
59          _RL tmpFac
60    #endif /* ALLOW_DIAGNOSTICS */
61    
62          IF ( usingPCoords ) THEN
63           ks        = Nr
64          ELSE
65           ks        = 1
66          ENDIF
67    
68    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
69    
70          IF ( doThetaClimRelax .OR. doSaltClimRelax ) THEN
71    C--   Start with surface restoring term :
72    
73           DO j = jMin, jMax
74            DO i = iMin, iMax
75    #ifdef ALLOW_SEAICE
76    C     Do not restore under sea-ice
77    C     Heat Flux (restoring term) :
78              surfaceForcingT(i,j,bi,bj) =
79         &      -lambdaThetaClimRelax(i,j,bi,bj) * (1-AREA(i,j,1,bi,bj))
80         &         *(theta(i,j,ks,bi,bj)-SST(i,j,bi,bj))
81         &         *drF(ks)*_hFacC(i,j,ks,bi,bj)
82    C     Salt Flux (restoring term) :
83              surfaceForcingS(i,j,bi,bj) =
84         &      -lambdaSaltClimRelax(i,j,bi,bj) * (1-AREA(i,j,1,bi,bj))
85         &         *(salt(i,j,ks,bi,bj)-SSS(i,j,bi,bj))
86         &         *drF(ks)*_hFacC(i,j,ks,bi,bj)
87    #else /* ifndef ALLOW_SEAICE */
88    C     Heat Flux (restoring term) :
89              surfaceForcingT(i,j,bi,bj) =
90         &      -lambdaThetaClimRelax(i,j,bi,bj)
91         &         *(theta(i,j,ks,bi,bj)-SST(i,j,bi,bj))
92         &         *drF(ks)*_hFacC(i,j,ks,bi,bj)
93    C     Salt Flux (restoring term) :
94              surfaceForcingS(i,j,bi,bj) =
95         &      -lambdaSaltClimRelax(i,j,bi,bj)
96         &         *(salt(i,j,ks,bi,bj)-SSS(i,j,bi,bj))
97         &         *drF(ks)*_hFacC(i,j,ks,bi,bj)
98    #endif /* ALLOW_SEAICE */
99            ENDDO
100           ENDDO
101    
102    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
103    #ifdef NONLIN_FRSURF
104    C-    T,S surface forcing will be applied (thermodynamics) after the update
105    C     of surf.thickness (hFac): account for change in surf.thickness
106           IF (staggerTimeStep.AND.nonlinFreeSurf.GT.0) THEN
107            IF (select_rStar.GT.0) THEN
108    # ifndef DISABLE_RSTAR_CODE
109             DO j=jMin,jMax
110              DO i=iMin,iMax
111                surfaceForcingT(i,j,bi,bj) = surfaceForcingT(i,j,bi,bj)
112         &                                  * rStarExpC(i,j,bi,bj)
113                surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
114         &                                  * rStarExpC(i,j,bi,bj)
115              ENDDO
116             ENDDO
117    # endif /* DISABLE_RSTAR_CODE */
118            ELSE
119             DO j=jMin,jMax
120              DO i=iMin,iMax
121               IF (ks.EQ.ksurfC(i,j,bi,bj)) THEN
122                surfaceForcingT(i,j,bi,bj) = surfaceForcingT(i,j,bi,bj)
123         &             *_recip_hFacC(i,j,ks,bi,bj)*hFac_surfC(i,j,bi,bj)
124                surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
125         &             *_recip_hFacC(i,j,ks,bi,bj)*hFac_surfC(i,j,bi,bj)
126               ENDIF
127              ENDDO
128             ENDDO
129            ENDIF
130           ENDIF
131    #endif /* NONLIN_FRSURF */
132    
133    #ifdef ALLOW_DIAGNOSTICS
134           IF ( useDiagnostics ) THEN
135    
136        if ( buoyancyRelation .eq. 'OCEANICP' ) then  C     tRelax (temperature relaxation [W/m2], positive <-> increasing Theta)
137         kSurface        = Nr          tmpFac = HeatCapacity_Cp*recip_horiVertRatio*rhoConst
138         convertVol2Mass = gravity*rhoConstFresh          CALL DIAGNOSTICS_SCALE_FILL(
139        else       &           surfaceForcingT(1-OLx,1-OLy,bi,bj),tmpFac,1,
140         kSurface        = 1       &                             'TRELAX  ',0, 1,2,bi,bj,myThid)
141         convertVol2Mass = 1. _d 0  
142        endif  C     sRelax (salt relaxation [g/m2/s], positive <-> increasing Salt)
143            tmpFac = recip_horiVertRatio*rhoConst
144            CALL DIAGNOSTICS_SCALE_FILL(
145         &           surfaceForcingS(1-OLx,1-OLy,bi,bj),tmpFac,1,
146         &                             'SRELAX  ',0, 1,2,bi,bj,myThid)
147    
148           ENDIF
149    #endif /* ALLOW_DIAGNOSTICS */
150    
151          ELSE
152    C--   No restoring for T & S : set surfaceForcingT,S to zero :
153    
154           DO j = jMin, jMax
155            DO i = iMin, iMax
156              surfaceForcingT(i,j,bi,bj) = 0. _d 0
157              surfaceForcingS(i,j,bi,bj) = 0. _d 0
158            ENDDO
159           ENDDO
160    
161    C--   end restoring / no restoring block.
162          ENDIF
163    
164    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
165    
166    C--   Surface Fluxes :
167    
168        DO j = jMin, jMax        DO j = jMin, jMax
169           DO i = iMin, iMax           DO i = iMin, iMax
170    
171  c     Zonal wind stress fu:  C     Zonal wind stress fu:
172            surfaceTendencyU(i,j,bi,bj) =            surfaceForcingU(i,j,bi,bj) =
173       &      fu(i,j,bi,bj)*horiVertRatio*recip_rhoConst       &      fu(i,j,bi,bj)*horiVertRatio*recip_rhoConst
174       &           *recip_drF(kSurface)*recip_hFacW(i,j,kSurface,bi,bj)  C     Meridional wind stress fv:
175  c     Meridional wind stress fv:            surfaceForcingV(i,j,bi,bj) =
           surfaceTendencyV(i,j,bi,bj) =  
176       &      fv(i,j,bi,bj)*horiVertRatio*recip_rhoConst       &      fv(i,j,bi,bj)*horiVertRatio*recip_rhoConst
177       &           *recip_drF(kSurface)*recip_hFacS(i,j,kSurface,bi,bj)  C     Net heat flux Qnet:
178  c     Net heat flux Qnet:            surfaceForcingT(i,j,bi,bj) = surfaceForcingT(i,j,bi,bj)
179            surfaceTendencyT(i,j,bi,bj) =       &       - ( Qnet(i,j,bi,bj)
180       &      -Qnet(i,j,bi,bj)*recip_Cp*horiVertRatio*recip_rhoConst  #ifdef SHORTWAVE_HEATING
181       &           *recip_drF(kSurface)*recip_hFacC(i,j,kSurface,bi,bj)       &          -Qsw(i,j,bi,bj)
182       &      -lambdaThetaClimRelax  #endif
183       &         *(theta(i,j,kSurface,bi,bj)-SST(i,j,bi,bj))       &         ) *recip_Cp*horiVertRatio*recip_rhoConst
184  C     Salt Flux (restoring term) :  C     Net Salt Flux :
185  C         surfaceTendencyS(i,j,bi,bj) =            surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
186  C    &      -lambdaSaltClimRelax*(salt(i,j,kSurface,bi,bj)-SSS(i,j,bi,bj))       &      -saltFlux(i,j,bi,bj)*horiVertRatio*recip_rhoConst
 C notes : because truncation is different when this tendency is splitted  
 C   in 2 parts, keep this salt flux with freshwater flux (see below)  
   
 #ifdef ALLOW_PASSIVE_TRACER  
 c ***  define the tracer surface tendency here ***  
 #endif /* ALLOW_PASSIVE_TRACER */  
187    
188           ENDDO           ENDDO
189        ENDDO        ENDDO
190    
191  c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
192  C     Surface salinity tendency and freshwater flux EmPmR:  C--   Fresh-water flux
193    
194    C-    Apply mask on Fresh-Water flux
195    C      (needed for SSH forcing, whether or not exactConserv is used)
196          IF ( useRealFreshWaterFlux ) THEN
197            DO j=1-OLy,sNy+OLy
198             DO i=1-OLx,sNx+OLx
199               EmPmR(i,j,bi,bj) = EmPmR(i,j,bi,bj)*maskH(i,j,bi,bj)
200             ENDDO
201            ENDDO
202          ENDIF
203    
204    #ifdef EXACT_CONSERV
205    c NB: synchronous time step: PmEpR lag 1 time step behind EmPmR
206    c     to stay consitent with volume change (=d/dt etaH).
207          IF ( staggerTimeStep ) THEN
208            DO j=1-OLy,sNy+OLy
209             DO i=1-OLx,sNx+OLx
210               PmEpR(i,j,bi,bj) = -EmPmR(i,j,bi,bj)
211             ENDDO
212            ENDDO
213          ENDIF
214    
215          IF ( (nonlinFreeSurf.GT.0 .OR. usingPCoords)
216         &     .AND. useRealFreshWaterFlux ) THEN
217    
218    c-  NonLin_FrSurf and RealFreshWaterFlux : PmEpR effectively changes
219    c   the water column height ; temp., salt, (tracer) flux associated
220    c   with this input/output of water is added here to the surface tendency.
221    c
222    
223           IF (temp_EvPrRn.NE.UNSET_RL) THEN
224            DO j = jMin, jMax
225             DO i = iMin, iMax
226              surfaceForcingT(i,j,bi,bj) = surfaceForcingT(i,j,bi,bj)
227         &      + PmEpR(i,j,bi,bj)
228         &         *( temp_EvPrRn - theta(i,j,ks,bi,bj) )
229         &         *convertEmP2rUnit
230             ENDDO
231            ENDDO
232           ENDIF
233    
234           IF (salt_EvPrRn.NE.UNSET_RL) THEN
235            DO j = jMin, jMax
236             DO i = iMin, iMax
237              surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
238         &      + PmEpR(i,j,bi,bj)
239         &         *( salt_EvPrRn - salt(i,j,ks,bi,bj) )
240         &         *convertEmP2rUnit
241             ENDDO
242            ENDDO
243           ENDIF
244    
245        IF (.NOT.useRealFreshWaterFlux .OR. nonlinFreeSurf .LE. 0 ) THEN  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
246          ELSE
247    #else /* EXACT_CONSERV */
248          IF (.TRUE.) THEN
249    #endif /* EXACT_CONSERV */
250    
251  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)
252  c   and is converted to a salt tendency.  c   and is converted to a salt tendency.
# Line 96  c   and is converted to a salt tendency. Line 255  c   and is converted to a salt tendency.
255  c- converts EmPmR to salinity tendency using surface local salinity  c- converts EmPmR to salinity tendency using surface local salinity
256          DO j = jMin, jMax          DO j = jMin, jMax
257           DO i = iMin, iMax           DO i = iMin, iMax
258            surfaceTendencyS(i,j,bi,bj) =            surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
259       &      + EmPmR(i,j,bi,bj)*salt(i,j,kSurface,bi,bj)       &      + EmPmR(i,j,bi,bj)*salt(i,j,ks,bi,bj)
260       &           *recip_drF(kSurface)*recip_hFacC(i,j,kSurface,bi,bj)       &         *convertEmP2rUnit
      &         *convertVol2Mass  
      &      -lambdaSaltClimRelax  
      &         *(salt(i,j,kSurface,bi,bj)-SSS(i,j,bi,bj))  
261           ENDDO           ENDDO
262          ENDDO          ENDDO
263         ELSE         ELSE
264  c- converts EmPmR to virtual salt flux using uniform salinity (default=35)  c- converts EmPmR to virtual salt flux using uniform salinity (default=35)
265          DO j = jMin, jMax          DO j = jMin, jMax
266           DO i = iMin, iMax           DO i = iMin, iMax
267            surfaceTendencyS(i,j,bi,bj) =            surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
268       &      + EmPmR(i,j,bi,bj)*convertFW2Salt       &      + EmPmR(i,j,bi,bj)*convertFW2Salt
269       &           *recip_drF(kSurface)*recip_hFacC(i,j,kSurface,bi,bj)       &         *convertEmP2rUnit
      &         *convertVol2Mass  
      &      -lambdaSaltClimRelax  
      &         *(salt(i,j,kSurface,bi,bj)-SSS(i,j,bi,bj))  
270           ENDDO           ENDDO
271          ENDDO          ENDDO
272         ENDIF         ENDIF
273    
274  #ifdef NONLIN_FRSURF        ENDIF
 c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  
       ELSE  
275    
276  c     Salt Flux (restoring term) :  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
        DO j = jMin, jMax  
          DO i = iMin, iMax  
           surfaceTendencyS(i,j,bi,bj) =  
      &      -lambdaSaltClimRelax  
      &         *(salt(i,j,kSurface,bi,bj)-SSS(i,j,bi,bj))  
          ENDDO  
        ENDDO  
277    
278  c-  NonLin_FrSurf and RealFreshWaterFlux : PmEpR effectively changes  #ifdef ALLOW_PTRACERS
279  c   the water column height ; temp., salt, (tracer) flux associated        IF ( usePTRACERS ) THEN
280  c   with this input/output of water is added here to the surface tendency.           CALL PTRACERS_FORCING_SURF(
281  c       I        bi, bj, iMin, iMax, jMin, jMax,
282  c NB: PmEpR lag 1 time step behind EmPmR ( PmEpR_n = - EmPmR_n-1 ) to stay       I        myTime,myIter,myThid )
283  c     consitent with volume change (=d/dt etaN).        ENDIF
284    #endif /* ALLOW_PTRACERS */
285    
286         IF (temp_EvPrRn.NE.UNSET_RL) THEN  #ifdef ATMOSPHERIC_LOADING
287    
288    C-- Atmospheric surface Pressure loading :
289    
290          IF ( usingZCoords ) THEN
291           IF ( useRealFreshWaterFlux ) THEN
292          DO j = jMin, jMax          DO j = jMin, jMax
293           DO i = iMin, iMax           DO i = iMin, iMax
294            surfaceTendencyT(i,j,bi,bj) = surfaceTendencyT(i,j,bi,bj)            phi0surf(i,j,bi,bj) = ( pload(i,j,bi,bj)
295       &      + PmEpR(i,j,bi,bj)       &                        +sIceLoad(i,j,bi,bj)*gravity
296       &         *( temp_EvPrRn - theta(i,j,kSurface,bi,bj) )       &                          )*recip_rhoConst
      &           *recip_drF(kSurface)*recip_hFacC(i,j,kSurface,bi,bj)  
      &         *convertVol2Mass  
297           ENDDO           ENDDO
298          ENDDO          ENDDO
299         ENDIF         ELSE
   
        IF (salt_EvPrRn.NE.UNSET_RL) THEN  
300          DO j = jMin, jMax          DO j = jMin, jMax
301           DO i = iMin, iMax           DO i = iMin, iMax
302            surfaceTendencyS(i,j,bi,bj) = surfaceTendencyS(i,j,bi,bj)            phi0surf(i,j,bi,bj) = pload(i,j,bi,bj)*recip_rhoConst
      &      + PmEpR(i,j,bi,bj)  
      &         *( salt_EvPrRn - salt(i,j,kSurface,bi,bj) )  
      &           *recip_drF(kSurface)*recip_hFacC(i,j,kSurface,bi,bj)  
      &         *convertVol2Mass  
303           ENDDO           ENDDO
304          ENDDO          ENDDO
305         ENDIF         ENDIF
306          ELSEIF ( usingPCoords ) THEN
307    C-- This is a hack used to read phi0surf from a file (ploadFile)
308    C   instead of computing it from bathymetry & density ref. profile.
309    C   The true atmospheric P-loading is not yet implemented for P-coord
310    C   (requires time varying dP(Nr) like dP(k-bottom) with NonLin FS).
311            DO j = jMin, jMax
312             DO i = iMin, iMax
313              phi0surf(i,j,bi,bj) = pload(i,j,bi,bj)
314             ENDDO
315            ENDDO
316          ENDIF
317    
318  #ifdef ALLOW_PASSIVE_TRACER  #endif /* ATMOSPHERIC_LOADING */
 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 */  
319    
320  #endif /* NONLIN_FRSURF */  #ifdef ALLOW_SHELFICE
321          IF ( usingZCoords ) THEN
322           IF ( useSHELFICE) THEN
323            DO j = jMin, jMax
324             DO i = iMin, iMax
325              phi0surf(i,j,bi,bj) = phi0surf(i,j,bi,bj)
326         &         + shelficeLoadAnomaly(i,j,bi,bj)*recip_rhoConst
327             ENDDO
328            ENDDO
329           ENDIF
330        ENDIF        ENDIF
331    #endif /* ALLOW_SHELFICE */
332    
333    #ifdef ALLOW_EBM
334    c--    Values for surfaceForcingT, surfaceForcingS
335    c      are overwritten by those produced by EBM
336    cph    AD recomputation problems if these IF useEBM are used
337    cph      IF ( useEBM ) THEN
338           CALL EBM_FORCING_SURF(
339         I        bi, bj, iMin, iMax, jMin, jMax,
340         I        myTime,myIter,myThid )
341    cph      ENDIF
342    #endif
343    
344        RETURN        RETURN
345        END        END
Legend:
Removed from v.1.9  
changed lines
  Added in v.1.36
  ViewVC Help
Powered by ViewVC 1.1.22