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

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

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

revision 1.35 by stephd, Mon Dec 19 19:09:35 2005 UTC revision 1.47 by jmc, Thu Aug 23 19:12:10 2007 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"
 #ifdef ALLOW_EXF  
 # include "EXF_OPTIONS.h"  
 #endif  
6    
7  CBOP  CBOP
8  C     !ROUTINE: EXTERNAL_FORCING_U  C     !ROUTINE: EXTERNAL_FORCING_U
# Line 74  C--   Forcing term Line 71  C--   Forcing term
71       &                      myTime, myThid )       &                      myTime, myThid )
72  #endif /* ALLOW_FIZHI */  #endif /* ALLOW_FIZHI */
73    
74    #ifdef ALLOW_MYPACKAGE
75          IF ( useMYPACKAGE ) CALL MYPACKAGE_TENDENCY_APPLY_U(
76         &                      iMin,iMax, jMin,jMax, bi,bj, kLev,
77         &                      myTime, myThid )
78    #endif /* ALLOW_MYPACKAGE */
79    
80  C     Add windstress momentum impulse into the top-layer  C     Add windstress momentum impulse into the top-layer
81        IF ( kLev .EQ. kSurface ) THEN        IF ( kLev .EQ. kSurface ) THEN
82  c      DO j=1,sNy  c      DO j=1,sNy
# Line 82  C-jmc: Without CD-scheme, this is OK ; b Line 85  C-jmc: Without CD-scheme, this is OK ; b
85          DO i=1,sNx+1          DO i=1,sNx+1
86           gU(i,j,kLev,bi,bj) = gU(i,j,kLev,bi,bj)           gU(i,j,kLev,bi,bj) = gU(i,j,kLev,bi,bj)
87       &   +foFacMom*surfaceForcingU(i,j,bi,bj)       &   +foFacMom*surfaceForcingU(i,j,bi,bj)
88       &   *recip_drF(kLev)*recip_hFacW(i,j,kLev,bi,bj)       &   *recip_drF(kLev)*_recip_hFacW(i,j,kLev,bi,bj)
89          ENDDO          ENDDO
90         ENDDO         ENDDO
91        ENDIF        ENDIF
# Line 172  C--   Forcing term Line 175  C--   Forcing term
175       &                      myTime, myThid )       &                      myTime, myThid )
176  #endif /* ALLOW_FIZHI */  #endif /* ALLOW_FIZHI */
177    
178    #ifdef ALLOW_MYPACKAGE
179          IF ( useMYPACKAGE ) CALL MYPACKAGE_TENDENCY_APPLY_V(
180         &                      iMin,iMax, jMin,jMax, bi,bj, kLev,
181         &                      myTime, myThid )
182    #endif /* ALLOW_MYPACKAGE */
183    
184  C     Add windstress momentum impulse into the top-layer  C     Add windstress momentum impulse into the top-layer
185        IF ( kLev .EQ. kSurface ) THEN        IF ( kLev .EQ. kSurface ) THEN
186         DO j=1,sNy+1         DO j=1,sNy+1
# Line 180  C-jmc: Without CD-scheme, this is OK ; b Line 189  C-jmc: Without CD-scheme, this is OK ; b
189          DO i=0,sNx+1          DO i=0,sNx+1
190           gV(i,j,kLev,bi,bj) = gV(i,j,kLev,bi,bj)           gV(i,j,kLev,bi,bj) = gV(i,j,kLev,bi,bj)
191       &   +foFacMom*surfaceForcingV(i,j,bi,bj)       &   +foFacMom*surfaceForcingV(i,j,bi,bj)
192       &   *recip_drF(kLev)*recip_hFacS(i,j,kLev,bi,bj)       &   *recip_drF(kLev)*_recip_hFacS(i,j,kLev,bi,bj)
193          ENDDO          ENDDO
194         ENDDO         ENDDO
195        ENDIF        ENDIF
# Line 228  C     == Global data == Line 237  C     == Global data ==
237  #include "GRID.h"  #include "GRID.h"
238  #include "DYNVARS.h"  #include "DYNVARS.h"
239  #include "FFIELDS.h"  #include "FFIELDS.h"
240    #include "SURFACE.h"
241    
242  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
243  C     == Routine arguments ==  C     == Routine arguments ==
# Line 277  C--   Forcing term Line 287  C--   Forcing term
287       &                      myTime, myThid )       &                      myTime, myThid )
288  #endif /* ALLOW_FIZHI */  #endif /* ALLOW_FIZHI */
289    
290    #ifdef ALLOW_MYPACKAGE
291          IF ( useMYPACKAGE ) CALL MYPACKAGE_TENDENCY_APPLY_T(
292         &                      iMin,iMax, jMin,jMax, bi,bj, kLev,
293         &                      myTime, myThid )
294    #endif /* ALLOW_MYPACKAGE */
295    
296  C     Add heat in top-layer  C     Add heat in top-layer
297        IF ( kLev .EQ. kSurface ) THEN        IF ( kLev .EQ. kSurface ) THEN
298         DO j=1,sNy         DO j=1,sNy
299          DO i=1,sNx          DO i=1,sNx
300           gT(i,j,kLev,bi,bj)=gT(i,j,kLev,bi,bj)           gT(i,j,kLev,bi,bj)=gT(i,j,kLev,bi,bj)
301       &     +surfaceForcingT(i,j,bi,bj)       &     +surfaceForcingT(i,j,bi,bj)
302       &     *recip_drF(kLev)*recip_hFacC(i,j,kLev,bi,bj)       &     *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
303          ENDDO          ENDDO
304         ENDDO         ENDDO
305        ENDIF        ENDIF
306    
307    cph#ifndef ALLOW_AUTODIFF_TAMC
308    cph I didnt put this ifndef here.
309          IF (linFSConserveTr) THEN
310           DO j=1,sNy
311            DO i=1,sNx
312              IF (kLev .EQ. ksurfC(i,j,bi,bj)) THEN
313                gT(i,j,kLev,bi,bj)=gT(i,j,kLev,bi,bj)
314         &        +TsurfCor*recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
315              ENDIF
316            ENDDO
317           ENDDO
318          ENDIF
319    cph#endif /* ndfef ALLOW_AUTODIFF_TAMC */
320    
321    #ifdef ALLOW_SHELFICE
322          IF ( useShelfIce )
323         &     CALL SHELFICE_FORCING_T(
324         I     iMin,iMax, jMin,jMax, bi,bj, kLev,
325         I     myTime, myThid )
326    #endif /* ALLOW_SHELFICE */
327    
328  #ifdef SHORTWAVE_HEATING  #ifdef SHORTWAVE_HEATING
329  C Penetrating SW radiation  C Penetrating SW radiation
330  c     IF ( usePenetratingSW ) THEN  c     IF ( usePenetratingSW ) THEN
331         swfracb(1)=abs(rF(klev))         swfracb(1)=abs(rF(klev))
332         swfracb(2)=abs(rF(klev+1))         swfracb(2)=abs(rF(klev+1))
333         CALL SWFRAC(         CALL SWFRAC(
334       I     two,minusone,       I             two, minusone,
335       I     myTime,myThid,       U             swfracb,
336       U     swfracb)       I             myTime, 1, myThid )
337         kp1 = klev+1         kp1 = klev+1
338         IF (klev.EQ.Nr) THEN         IF (klev.EQ.Nr) THEN
339          kp1 = klev          kp1 = klev
# Line 307  c     IF ( usePenetratingSW ) THEN Line 344  c     IF ( usePenetratingSW ) THEN
344           gT(i,j,klev,bi,bj) = gT(i,j,klev,bi,bj)           gT(i,j,klev,bi,bj) = gT(i,j,klev,bi,bj)
345       &   -Qsw(i,j,bi,bj)*(swfracb(1)*maskC(i,j,klev,bi,bj)       &   -Qsw(i,j,bi,bj)*(swfracb(1)*maskC(i,j,klev,bi,bj)
346       &                   -swfracb(2)*maskC(i,j,kp1, bi,bj))       &                   -swfracb(2)*maskC(i,j,kp1, bi,bj))
347       &    *recip_Cp*recip_rhoConst       &    *recip_Cp*mass2rUnit
348       &    *recip_drF(klev)*recip_hFacC(i,j,kLev,bi,bj)       &    *recip_drF(klev)*_recip_hFacC(i,j,kLev,bi,bj)
349          ENDDO          ENDDO
350         ENDDO         ENDDO
351  c     ENDIF  c     ENDIF
352  #endif  #endif
353    
354  #ifdef ALLOW_RBCS  #ifdef ALLOW_RBCS
355         if (useRBCS) then         IF (useRBCS) THEN
356            call RBCS_ADD_TENDENCY(bi,bj,klev, 1,            CALL RBCS_ADD_TENDENCY(bi,bj,klev, 1,
357       &                            myTime, myThid )       &                            myTime, myThid )
        endif  
 #endif  
   
 #ifdef ALLOW_CLIMTEMP_RELAXATION  
        IF ( tauThetaClimRelax3Dim .NE. 0. ) THEN  
         DO j=1,sNy  
          DO i=1,sNx  
           gT(i,j,klev,bi,bj) = gT(i,j,klev,bi,bj)  
      &     -1./tauThetaClimRelax3Dim  
      &         *(theta(i,j,klev,bi,bj)-thetaStar(i,j,klev,bi,bj))  
      &         *hFacC(i,j,klev,bi,bj)*recip_hFacC(i,j,kLev,bi,bj)  
          ENDDO  
         ENDDO  
358         ENDIF         ENDIF
359  #endif  #endif
360    
# Line 372  C     == Global data == Line 396  C     == Global data ==
396  #include "GRID.h"  #include "GRID.h"
397  #include "DYNVARS.h"  #include "DYNVARS.h"
398  #include "FFIELDS.h"  #include "FFIELDS.h"
399    #include "SURFACE.h"
400    #ifdef ALLOW_SALT_PLUME
401    #ifdef ALLOW_SEAICE
402    #include "SEAICE_PARAMS.h"
403    #endif /* ALLOW_SEAICE */
404    #endif /* ALLOW_SALT_PLUME */
405    
406  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
407  C     == Routine arguments ==  C     == Routine arguments ==
# Line 392  C     kSurface  :: index of surface laye Line 422  C     kSurface  :: index of surface laye
422        INTEGER i, j        INTEGER i, j
423        INTEGER kSurface        INTEGER kSurface
424  CEOP  CEOP
425    #ifdef ALLOW_SALT_PLUME
426          _RL saltPlume
427    #endif /* ALLOW_SALT_PLUME */
428    
429        IF ( fluidIsAir ) THEN        IF ( fluidIsAir ) THEN
430         kSurface = 0         kSurface = 0
# Line 414  C--   Forcing term Line 447  C--   Forcing term
447       &                      myTime, myThid )       &                      myTime, myThid )
448  #endif /* ALLOW_FIZHI */  #endif /* ALLOW_FIZHI */
449    
450    #ifdef ALLOW_MYPACKAGE
451          IF ( useMYPACKAGE ) CALL MYPACKAGE_TENDENCY_APPLY_S(
452         &                      iMin,iMax, jMin,jMax, bi,bj, kLev,
453         &                      myTime, myThid )
454    #endif /* ALLOW_MYPACKAGE */
455    
456  C     Add fresh-water in top-layer  C     Add fresh-water in top-layer
457        IF ( kLev .EQ. kSurface ) THEN        IF ( kLev .EQ. kSurface ) THEN
458         DO j=1,sNy         DO j=1,sNy
459          DO i=1,sNx          DO i=1,sNx
460           gS(i,j,kLev,bi,bj)=gS(i,j,kLev,bi,bj)           gS(i,j,kLev,bi,bj)=gS(i,j,kLev,bi,bj)
461       &     +surfaceForcingS(i,j,bi,bj)       &     +surfaceForcingS(i,j,bi,bj)
462       &     *recip_drF(kLev)*recip_hFacC(i,j,kLev,bi,bj)       &     *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
463          ENDDO          ENDDO
464         ENDDO         ENDDO
465        ENDIF        ENDIF
466    
467  #ifdef ALLOW_RBCS  cph#ifndef ALLOW_AUTODIFF_TAMC
468         if (useRBCS) then  cph I didnt put this ifndef here.
469            call RBCS_ADD_TENDENCY(bi,bj,klev, 2,        IF (linFSConserveTr) THEN
470       &                            myTime, myThid )         DO j=1,sNy
471         endif          DO i=1,sNx
472  #endif            IF (kLev .EQ. ksurfC(i,j,bi,bj)) THEN
473                gS(i,j,kLev,bi,bj)=gS(i,j,kLev,bi,bj)
474         &        +SsurfCor*recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
475              ENDIF
476            ENDDO
477           ENDDO
478          ENDIF
479    cph#endif /* ndfef ALLOW_AUTODIFF_TAMC */
480    
481  #ifdef ALLOW_CLIMSALT_RELAXATION  #ifdef ALLOW_SHELFICE
482         IF ( tauSaltClimRelax3Dim .NE. 0. ) THEN        IF ( useShelfIce )
483          DO j=1,sNy       &     CALL SHELFICE_FORCING_S(
484           DO i=1,sNx       I     iMin,iMax, jMin,jMax, bi,bj, kLev,
485            gS(i,j,klev,bi,bj) = gS(i,j,klev,bi,bj)       I     myTime, myThid )
486       &     -1./tauSaltClimRelax3Dim  #endif /* ALLOW_SHELFICE */
487       &         *(salt(i,j,klev,bi,bj)-saltStar(i,j,klev,bi,bj))  
488       &         *hFacC(i,j,klev,bi,bj)*recip_hFacC(i,j,kLev,bi,bj)  #ifdef ALLOW_SALT_PLUME
489           ENDDO  C saltPlume is the amount of salt rejected by ice while freezing;
490    C it is here redistributed to multiple vertical levels as per
491    C Duffy et al. (GRL 1999)
492           DO j=1,sNy
493            DO i=1,sNx
494              saltPlume = 0.
495    #ifdef ALLOW_SEAICE
496              IF ( saltFlux(i,j,bi,bj) .GT. 0. .AND.
497         &         salt(i,j,kSurface,bi,bj)  .GT. SEAICE_salinity ) THEN
498               saltPlume = (salt(i,j,kSurface,bi,bj)-SEAICE_salinity) *
499         &          saltFlux(i,j,bi,bj) / salt(i,j,kSurface,bi,bj)
500              ENDIF
501    #endif /* ALLOW_SEAICE */
502              IF ( SaltPlumeDepth(i,j,bi,bj) .GT. -rF(kLev) ) THEN
503               gS(i,j,kLev,bi,bj)=gS(i,j,kLev,bi,bj)
504         &          +saltPlume*mass2rUnit
505         &          *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
506         &          *min(drF(kLev),SaltPlumeDepth(i,j,bi,bj)+rF(kLev))
507         &          /SaltPlumeDepth(i,j,bi,bj)
508              ENDIF
509          ENDDO          ENDDO
510           ENDDO
511    #endif /* ALLOW_SALT_PLUME */
512    
513    #ifdef ALLOW_RBCS
514           IF (useRBCS) THEN
515              CALL RBCS_ADD_TENDENCY(bi,bj,klev, 2,
516         &                            myTime, myThid )
517         ENDIF         ENDIF
518  #endif  #endif /* ALLOW_RBCS */
519    
520  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
521        IF (useOBCS) THEN        IF (useOBCS) THEN
# Line 451  C     Add fresh-water in top-layer Line 523  C     Add fresh-water in top-layer
523       I           iMin,iMax, jMin,jMax, bi,bj, kLev,       I           iMin,iMax, jMin,jMax, bi,bj, kLev,
524       I           myTime, myThid )       I           myTime, myThid )
525        ENDIF        ENDIF
526  #endif  #endif /* ALLOW_OBCS */
527    
528        RETURN        RETURN
529        END        END
Legend:
Removed from v.1.35  
changed lines
  Added in v.1.47
  ViewVC Help
Powered by ViewVC 1.1.22