/[MITgcm]/MITgcm/verification/hs94.cs-32x32x5/code/external_forcing.F
ViewVC logotype

Diff of /MITgcm/verification/hs94.cs-32x32x5/code/external_forcing.F

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

revision 1.1.2.1 by adcroft, Mon Apr 9 20:01:16 2001 UTC revision 1.8 by jmc, Wed Aug 20 20:23:10 2014 UTC
# Line 3  C $Name$ Line 3  C $Name$
3    
4  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
5    
6  CStartOfInterface  CBOP
7    C     !ROUTINE: EXTERNAL_FORCING_U
8    C     !INTERFACE:
9        SUBROUTINE EXTERNAL_FORCING_U(        SUBROUTINE EXTERNAL_FORCING_U(
10       I           iMin, iMax, jMin, jMax,bi,bj,kLev,       I           iMin,iMax, jMin,jMax, bi,bj, kLev,
11       I           myCurrentTime,myThid)       I           myTime, myThid )
12  C     /==========================================================\  C     !DESCRIPTION: \bv
13  C     | S/R EXTERNAL_FORCING_U                                   |  C     *==========================================================*
14  C     | o Contains problem specific forcing for zonal velocity.  |  C     | S/R EXTERNAL_FORCING_U
15  C     |==========================================================|  C     | o Contains problem specific forcing for zonal velocity.
16  C     | Adds terms to gU for forcing by external sources         |  C     *==========================================================*
17  C     | e.g. wind stress, bottom friction etc..................  |  C     | Adds terms to gU for forcing by external sources
18  C     \==========================================================/  C     | e.g. wind stress, bottom friction etc ...
19        IMPLICIT NONE  C     *==========================================================*
20    C     \ev
21    
22    C     !USES:
23          IMPLICIT NONE
24  C     == Global data ==  C     == Global data ==
25  #include "SIZE.h"  #include "SIZE.h"
26  #include "EEPARAMS.h"  #include "EEPARAMS.h"
27  #include "PARAMS.h"  #include "PARAMS.h"
28  #include "GRID.h"  #include "GRID.h"
29    #include "SURFACE.h"
30  #include "DYNVARS.h"  #include "DYNVARS.h"
31  #include "FFIELDS.h"  #include "FFIELDS.h"
32    
33    C     !INPUT/OUTPUT PARAMETERS:
34  C     == Routine arguments ==  C     == Routine arguments ==
35  C     iMin - Working range of tile for applying forcing.  C     iMin,iMax :: Working range of x-index for applying forcing.
36  C     iMax  C     jMin,jMax :: Working range of y-index for applying forcing.
37  C     jMin  C     bi,bj     :: Current tile indices
38  C     jMax  C     kLev      :: Current vertical level index
39  C     kLev  C     myTime    :: Current time in simulation
40    C     myThid    :: Thread Id number
41        INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj        INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
42        _RL myCurrentTime        _RL myTime
43        INTEGER myThid        INTEGER myThid
 CEndOfInterface  
44    
45    C     !LOCAL VARIABLES:
46  C     == Local variables ==  C     == Local variables ==
47  C     Loop counters  C     i,j       :: Loop counters
48        INTEGER I, J        INTEGER i, j
49  C     _RL uKf  CEOP
50  C     _RL levelOfGround        _RL recip_P0g, termP, rFullDepth
51  C     _RL criticalLevel        _RL kV, kF, sigma_b
 C     _RL levelOfVelPoint  
 C     _RL dist1  
 C     _RL dist2  
 C     _RL decayFac  
 C     _RL velDragHeightFac  
       _RL termP,kV,kF  
52    
53  C--   Forcing term(s)  C--   Forcing term(s)
54        kF=1./86400.        kF=1. _d 0/86400. _d 0
55        DO J=jMin,jMax        sigma_b = 0.7 _d 0
56         DO I=iMin,iMax        rFullDepth = rF(1)-rF(Nr+1)
57          IF ( HFacW(i,j,kLev,bi,bj) .GT. 0. ) THEN  c     DO j=1,sNy
58  C        termP=0.5*( rF(kLev) + min( rF(kLev+1) ,  C-jmc: Without CD-scheme, this is OK ; but with CD-scheme, needs to cover [0:sNy+1]
59  C    &           min(H(I,J,bi,bj),H(I,J-1,bi,bj))            ) )        DO j=0,sNy+1
60           termP=0.5*( rF(kLev) + rF(kLev+1) )         DO i=1,sNx+1
61           kV=kF*MAX(0., (termP*recip_Rcol(I,J,bi,bj)-0.7)/(1.-0.7) )          IF ( hFacW(i,j,kLev,bi,bj) .GT. 0. ) THEN
62             IF ( selectSigmaCoord.EQ.0 ) THEN
63              recip_P0g=MAX(recip_Rcol(i,j,bi,bj),recip_Rcol(i-1,j,bi,bj))
64              termP=0.5 _d 0*( MIN(rF(kLev)*recip_P0g,1. _d 0)
65         &                    +rF(kLev+1)*recip_P0g )
66             ELSE
67    C-- Pressure at U.point :
68    c         midP = rLowW(i,j,bi,bj) + aHybSigmC(k)*rFullDepth
69    c    &         + bHybSigmC(k)
70    c    &          *(etaHw(i,j,bi,bj)+rSurfW(i,j,bi,bj)-rLowW(i,j,bi,bj))
71    C-- Sigma at U.point :
72    c         termP = ( midP - rLowW(i,j,bi,bj))
73    c    &          /(etaHw(i,j,bi,bj)+rSurfW(i,j,bi,bj)-rLowW(i,j,bi,bj))
74    C-  which simplifies to:
75              termP = aHybSigmC(kLev)*rFullDepth
76         &          /(etaHw(i,j,bi,bj)+rSurfW(i,j,bi,bj)-rLowW(i,j,bi,bj))
77         &          + bHybSigmC(kLev)
78             ENDIF
79             kV=kF*MAX( 0. _d 0, (termP-sigma_b)/(1. _d 0-sigma_b) )
80           gU(i,j,kLev,bi,bj)=gU(i,j,kLev,bi,bj)           gU(i,j,kLev,bi,bj)=gU(i,j,kLev,bi,bj)
81       &                      -kV*uVel(i,j,kLev,bi,bj)       &                     -kV*uVel(i,j,kLev,bi,bj)
82          ENDIF          ENDIF
83         ENDDO         ENDDO
84        ENDDO        ENDDO
85    
86        RETURN        RETURN
87        END        END
88  CStartOfInterface  
89    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
90    CBOP
91    C     !ROUTINE: EXTERNAL_FORCING_V
92    C     !INTERFACE:
93        SUBROUTINE EXTERNAL_FORCING_V(        SUBROUTINE EXTERNAL_FORCING_V(
94       I           iMin, iMax, jMin, jMax,bi,bj,kLev,       I           iMin,iMax, jMin,jMax, bi,bj, kLev,
95       I           myCurrentTime,myThid)       I           myTime, myThid )
96  C     /==========================================================\  C     !DESCRIPTION: \bv
97  C     | S/R EXTERNAL_FORCING_V                                   |  C     *==========================================================*
98  C     | o Contains problem specific forcing for merid velocity.  |  C     | S/R EXTERNAL_FORCING_V
99  C     |==========================================================|  C     | o Contains problem specific forcing for merid velocity.
100  C     | Adds terms to gV for forcing by external sources         |  C     *==========================================================*
101  C     | e.g. wind stress, bottom friction etc..................  |  C     | Adds terms to gV for forcing by external sources
102  C     \==========================================================/  C     | e.g. wind stress, bottom friction etc ...
103        IMPLICIT NONE  C     *==========================================================*
104    C     \ev
105    
106    C     !USES:
107          IMPLICIT NONE
108  C     == Global data ==  C     == Global data ==
109  #include "SIZE.h"  #include "SIZE.h"
110  #include "EEPARAMS.h"  #include "EEPARAMS.h"
111  #include "PARAMS.h"  #include "PARAMS.h"
112  #include "GRID.h"  #include "GRID.h"
113    #include "SURFACE.h"
114  #include "DYNVARS.h"  #include "DYNVARS.h"
115  #include "FFIELDS.h"  #include "FFIELDS.h"
116    
117    C     !INPUT/OUTPUT PARAMETERS:
118  C     == Routine arguments ==  C     == Routine arguments ==
119  C     iMin - Working range of tile for applying forcing.  C     iMin,iMax :: Working range of x-index for applying forcing.
120  C     iMax  C     jMin,jMax :: Working range of y-index for applying forcing.
121  C     jMin  C     bi,bj     :: Current tile indices
122  C     jMax  C     kLev      :: Current vertical level index
123  C     kLev  C     myTime    :: Current time in simulation
124    C     myThid    :: Thread Id number
125        INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj        INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
126        _RL myCurrentTime        _RL myTime
127        INTEGER myThid        INTEGER myThid
128  CEndOfInterface  
129    C     !LOCAL VARIABLES:
130  C     == Local variables ==  C     == Local variables ==
131  C     Loop counters  C     i,j       :: Loop counters
132        INTEGER I, J        INTEGER i, j
133  C     _RL uKf  CEOP
134  C     _RL levelOfGround        _RL recip_P0g, termP, rFullDepth
135  C     _RL criticalLevel        _RL kV, kF, sigma_b
 C     _RL levelOfVelPoint  
 C     _RL dist1  
 C     _RL dist2  
 C     _RL decayFac  
 C     _RL velDragHeightFac  
       _RL termP,kV,kF  
136    
137  C--   Forcing term(s)  C--   Forcing term(s)
138        kF=1./86400.        kF=1. _d 0/86400. _d 0
139        DO J=jMin,jMax        sigma_b = 0.7 _d 0
140         DO I=iMin,iMax        rFullDepth = rF(1)-rF(Nr+1)
141          IF ( HFacS(i,j,kLev,bi,bj) .GT. 0. ) THEN        DO j=1,sNy+1
142  C        termP=0.5*( rF(kLev) + min( rF(kLev+1) ,  c      DO i=1,sNx
143  C    &           min(H(I,J,bi,bj),H(I,J-1,bi,bj))            ) )  C-jmc: Without CD-scheme, this is OK ; but with CD-scheme, needs to cover [0:sNx+1]
144           termP=0.5*( rF(kLev) + rF(kLev+1) )         DO i=0,sNx+1
145           kV=kF*MAX(0., (termP*recip_Rcol(I,J,bi,bj)-0.7)/(1.-0.7) )          IF ( hFacS(i,j,kLev,bi,bj) .GT. 0. ) THEN
146             IF ( selectSigmaCoord.EQ.0 ) THEN
147              recip_P0g=MAX(recip_Rcol(i,j,bi,bj),recip_Rcol(i,j-1,bi,bj))
148              termP=0.5 _d 0*( MIN(rF(kLev)*recip_P0g,1. _d 0)
149         &                    +rF(kLev+1)*recip_P0g )
150             ELSE
151    C-- Pressure at V.point :
152    c         midP = rLowS(i,j,bi,bj) + aHybSigmC(k)*rFullDepth
153    c    &         + bHybSigmC(k)
154    c    &          *(etaHs(i,j,bi,bj)+rSurfS(i,j,bi,bj)-rLowS(i,j,bi,bj))
155    C-- Sigma at V.point :
156    c         termP = ( midP - rLowS(i,j,bi,bj))
157    c    &          /(etaHs(i,j,bi,bj)+rSurfS(i,j,bi,bj)-rLowS(i,j,bi,bj))
158    C-  which simplifies to:
159              termP = aHybSigmC(kLev)*rFullDepth
160         &          /(etaHs(i,j,bi,bj)+rSurfS(i,j,bi,bj)-rLowS(i,j,bi,bj))
161         &          + bHybSigmC(kLev)
162             ENDIF
163             kV=kF*MAX( 0. _d 0, (termP-sigma_b)/(1. _d 0-sigma_b) )
164           gV(i,j,kLev,bi,bj)=gV(i,j,kLev,bi,bj)           gV(i,j,kLev,bi,bj)=gV(i,j,kLev,bi,bj)
165       &                      -kV*vVel(i,j,kLev,bi,bj)       &                      -kV*vVel(i,j,kLev,bi,bj)
166          ENDIF          ENDIF
# Line 127  C    &           min(H(I,J,bi,bj),H(I,J- Line 169  C    &           min(H(I,J,bi,bj),H(I,J-
169    
170        RETURN        RETURN
171        END        END
172  CStartOfInterface  
173    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
174    CBOP
175    C     !ROUTINE: EXTERNAL_FORCING_T
176    C     !INTERFACE:
177        SUBROUTINE EXTERNAL_FORCING_T(        SUBROUTINE EXTERNAL_FORCING_T(
178       I           iMin, iMax, jMin, jMax,bi,bj,kLev,       I           iMin,iMax, jMin,jMax, bi,bj, kLev,
179       I           myCurrentTime,myThid)       I           myTime, myThid )
180  C     /==========================================================\  C     !DESCRIPTION: \bv
181  C     | S/R EXTERNAL_FORCING_T                                   |  C     *==========================================================*
182  C     | o Contains problem specific forcing for temperature.     |  C     | S/R EXTERNAL_FORCING_T
183  C     |==========================================================|  C     | o Contains problem specific forcing for temperature.
184  C     | Adds terms to gT for forcing by external sources         |  C     *==========================================================*
185  C     | e.g. heat flux, climatalogical relaxation..............  |  C     | Adds terms to gT for forcing by external sources
186  C     \==========================================================/  C     | e.g. heat flux, climatalogical relaxation, etc ...
187        IMPLICIT NONE  C     *==========================================================*
188    C     \ev
189    
190    C     !USES:
191          IMPLICIT NONE
192  C     == Global data ==  C     == Global data ==
193  #include "SIZE.h"  #include "SIZE.h"
194  #include "EEPARAMS.h"  #include "EEPARAMS.h"
# Line 148  C     == Global data == Line 197  C     == Global data ==
197  #include "DYNVARS.h"  #include "DYNVARS.h"
198  #include "FFIELDS.h"  #include "FFIELDS.h"
199    
200    C     !INPUT/OUTPUT PARAMETERS:
201  C     == Routine arguments ==  C     == Routine arguments ==
202  C     iMin - Working range of tile for applying forcing.  C     iMin,iMax :: Working range of x-index for applying forcing.
203  C     iMax  C     jMin,jMax :: Working range of y-index for applying forcing.
204  C     jMin  C     bi,bj     :: Current tile indices
205  C     jMax  C     kLev      :: Current vertical level index
206  C     kLev  C     myTime    :: Current time in simulation
207    C     myThid    :: Thread Id number
208        INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj        INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
209        _RL myCurrentTime        _RL myTime
210        INTEGER myThid        INTEGER myThid
 CEndOfInterface  
211    
212    C     !LOCAL VARIABLES:
213  C     == Local variables ==  C     == Local variables ==
214  C     Loop counters  C     i,j       :: Loop counters
215        INTEGER I, J        INTEGER i, j
216        _RL thetaLim,kT,ka,ks,term1,term2,thetaEq,termP,rSurf  CEOP
217          _RL thetaLim,kT,ka,ks,sigma_b,term1,term2,thetaEq
218          _RL termP, rFullDepth
219    
220  C--   Forcing term(s)  C--   Forcing term(s)
221        rSurf=1.E5        ka=1. _d 0/(40. _d 0*86400. _d 0)
222        ka=1./(40.*86400.)        ks=1. _d 0/(4. _d 0 *86400. _d 0)
223        ks=1./(4. *86400.)        sigma_b = 0.7 _d 0
224        DO J=jMin,jMax        rFullDepth = rF(1)-rF(Nr+1)
225         DO I=iMin,iMax        DO j=1,sNy
226         term1=60.*(sin(yC(I,J,bi,bj)*deg2rad)**2)         DO i=1,sNx
227  C      termP=0.5*( rF(kLev) + min( rF(kLev+1) , H(I,J,bi,bj) ) )           term1=60. _d 0*(SIN(yC(i,j,bi,bj)*deg2rad)**2)
228         termP=0.5*( rF(kLev) + rF(kLev+1) )           termP=0.5 _d 0*( rF(kLev) + rF(kLev+1) )
229         term2=10.*log(termP/rSurf)           term2=10. _d 0*LOG(termP/atm_po)
230       &          *(cos(yC(I,J,bi,bj)*deg2rad)**2)       &            *(COS(yC(i,j,bi,bj)*deg2rad)**2)
231         thetaLim = 200. / ((termP/rSurf)**(2./7.))           thetaLim = 200. _d 0/ ((termP/atm_po)**atm_kappa)
232         thetaEq=315.-term1-term2           thetaEq=315. _d 0-term1-term2
233         thetaEq=MAX(thetaLim,thetaEq)           thetaEq=MAX(thetaLim,thetaEq)
234          kT=ka+(ks-ka)           IF ( selectSigmaCoord.EQ.0 ) THEN
235       &    *MAX(0., (termP*recip_Rcol(I,J,bi,bj)-0.7)/(1.-0.7) )            termP=0.5 _d 0*( MIN(rF(kLev),Ro_surf(i,j,bi,bj))+rF(kLev+1) )
236       &    *COS((yC(I,J,bi,bj)*deg2rad))**4       &                  *recip_Rcol(i,j,bi,bj)
237             ELSE
238    C-- Pressure at T.point :
239    c         midP = R_low(i,j,bi,bj) + aHybSigmC(k)*rFullDepth
240    c    &         + bHybSigmC(k)
241    c    &          *(etaH(i,j,bi,bj)+Ro_surf(i,j,bi,bj)-R_low(i,j,bi,bj))
242    C-- Sigma at T.point :
243    c         termP = ( midP - R_low(i,j,bi,bj))
244    c    &          /(etaH(i,j,bi,bj)+Ro_surf(i,j,bi,bj)-R_low(i,j,bi,bj))
245    C-  which simplifies to:
246              termP = aHybSigmC(kLev)*rFullDepth
247         &          /(etaH(i,j,bi,bj)+Ro_surf(i,j,bi,bj)-R_low(i,j,bi,bj))
248         &          + bHybSigmC(kLev)
249             ENDIF
250             kT=ka+(ks-ka)
251         &     *MAX(0. _d 0, (termP-sigma_b)/(1. _d 0-sigma_b) )
252         &     *COS((yC(i,j,bi,bj)*deg2rad))**4
253           gT(i,j,kLev,bi,bj)=gT(i,j,kLev,bi,bj)           gT(i,j,kLev,bi,bj)=gT(i,j,kLev,bi,bj)
254       &        - kT*( theta(I,J,kLev,bi,bj)-thetaEq )       &        - kT*( theta(i,j,kLev,bi,bj)-thetaEq )
255       &            *maskC(i,j,kLev,bi,bj)       &            *maskC(i,j,kLev,bi,bj)
256         ENDDO         ENDDO
257        ENDDO        ENDDO
258    
259        RETURN        RETURN
260        END        END
261  CStartOfInterface  
262    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
263    CBOP
264    C     !ROUTINE: EXTERNAL_FORCING_S
265    C     !INTERFACE:
266        SUBROUTINE EXTERNAL_FORCING_S(        SUBROUTINE EXTERNAL_FORCING_S(
267       I           iMin, iMax, jMin, jMax,bi,bj,kLev,       I           iMin,iMax, jMin,jMax, bi,bj, kLev,
268       I           myCurrentTime,myThid)       I           myTime, myThid )
269  C     /==========================================================\  
270  C     | S/R EXTERNAL_FORCING_S                                   |  C     !DESCRIPTION: \bv
271  C     | o Contains problem specific forcing for merid velocity.  |  C     *==========================================================*
272  C     |==========================================================|  C     | S/R EXTERNAL_FORCING_S
273  C     | Adds terms to gS for forcing by external sources         |  C     | o Contains problem specific forcing for merid velocity.
274  C     | e.g. fresh-water flux, climatalogical relaxation.......  |  C     *==========================================================*
275  C     \==========================================================/  C     | Adds terms to gS for forcing by external sources
276        IMPLICIT NONE  C     | e.g. fresh-water flux, climatalogical relaxation, etc ...
277    C     *==========================================================*
278    C     \ev
279    
280    C     !USES:
281          IMPLICIT NONE
282  C     == Global data ==  C     == Global data ==
283  #include "SIZE.h"  #include "SIZE.h"
284  #include "EEPARAMS.h"  #include "EEPARAMS.h"
# Line 210  C     == Global data == Line 287  C     == Global data ==
287  #include "DYNVARS.h"  #include "DYNVARS.h"
288  #include "FFIELDS.h"  #include "FFIELDS.h"
289    
290    C     !INPUT/OUTPUT PARAMETERS:
291  C     == Routine arguments ==  C     == Routine arguments ==
292  C     iMin - Working range of tile for applying forcing.  C     iMin,iMax :: Working range of x-index for applying forcing.
293  C     iMax  C     jMin,jMax :: Working range of y-index for applying forcing.
294  C     jMin  C     bi,bj     :: Current tile indices
295  C     jMax  C     kLev      :: Current vertical level index
296  C     kLev  C     myTime    :: Current time in simulation
297    C     myThid    :: Thread Id number
298        INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj        INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
299        _RL myCurrentTime        _RL myTime
300        INTEGER myThid        INTEGER myThid
 CEndOfInterface  
301    
302    C     !LOCAL VARIABLES:
303  C     == Local variables ==  C     == Local variables ==
304  C     Loop counters  C     i,j       :: Loop counters
305        INTEGER I, J  c     INTEGER i, j
306    CEOP
307    
308  C--   Forcing term(s)  C--   Forcing term(s)
309    
Legend:
Removed from v.1.1.2.1  
changed lines
  Added in v.1.8
  ViewVC Help
Powered by ViewVC 1.1.22