/[MITgcm]/MITgcm/verification/hs94.128x64x5/code/external_forcing.F
ViewVC logotype

Diff of /MITgcm/verification/hs94.128x64x5/code/external_forcing.F

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

revision 1.2 by adcroft, Fri Feb 2 21:36:33 2001 UTC revision 1.8 by jmc, Sun Jul 17 18:40:10 2005 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    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"
# Line 23  C     == Global data == Line 29  C     == Global data ==
29  #include "DYNVARS.h"  #include "DYNVARS.h"
30  #include "FFIELDS.h"  #include "FFIELDS.h"
31    
32    C     !INPUT/OUTPUT PARAMETERS:
33  C     == Routine arguments ==  C     == Routine arguments ==
34  C     iMin - Working range of tile for applying forcing.  C     iMin,iMax :: Working range of x-index for applying forcing.
35  C     iMax  C     jMin,jMax :: Working range of y-index for applying forcing.
36  C     jMin  C     bi,bj     :: Current tile indices
37  C     jMax  C     kLev      :: Current vertical level index
38  C     kLev  C     myTime    :: Current time in simulation
39    C     myThid    :: Thread Id number
40        INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj        INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
41        _RL myCurrentTime        _RL myTime
42        INTEGER myThid        INTEGER myThid
 CEndOfInterface  
43    
44    C     !LOCAL VARIABLES:
45  C     == Local variables ==  C     == Local variables ==
46  C     Loop counters  C     i,j       :: Loop counters
47        INTEGER I, J        INTEGER i, j
48  C     _RL uKf  CEOP
49  C     _RL levelOfGround        _RL recip_P0g, termP, kV, kF, sigma_b
50  C     _RL criticalLevel  
51  C     _RL levelOfVelPoint  C--   Forcing term(s)
52  C     _RL dist1        kF=1. _d 0/86400. _d 0
53  C     _RL dist2        sigma_b = 0.7 _d 0
54  C     _RL decayFac  c     DO j=1,sNy
55  C     _RL velDragHeightFac  C-jmc: Without CD-scheme, this is OK ; but with CD-scheme, needs to cover [0:sNy+1]
56        _RL termP,kV,kF        DO j=0,sNy+1
57           DO i=1,sNx+1
58        kF=1./86400.          IF ( hFacW(i,j,kLev,bi,bj) .GT. 0. ) THEN
59        DO J=jMin,jMax           recip_P0g=MAX(recip_Rcol(i,j,bi,bj),recip_Rcol(i-1,j,bi,bj))
60         DO I=iMin,iMax           termP=0.5 _d 0*( MIN(rF(kLev)*recip_P0g,1. _d 0)
61          IF ( HFacW(i,j,kLev,bi,bj) .GT. 0. ) THEN       &                   +rF(kLev+1)*recip_P0g )
62  C        termP=0.5*( rF(kLev) + min( rF(kLev+1) ,  c        termP=0.5 _d 0*( rF(kLev) + rF(kLev+1) )*recip_P0g
63  C    &           min(H(I,J,bi,bj),H(I,J-1,bi,bj))            ) )           kV=kF*MAX( 0. _d 0, (termP-sigma_b)/(1. _d 0-sigma_b) )
          termP=0.5*( rF(kLev) + rF(kLev+1) )  
 C        termP=rC(kLev)  
          kV=kF*MAX(0., (termP*recip_H(I,J,bi,bj)-0.7)/(1.-0.7) )  
64           gU(i,j,kLev,bi,bj)=gU(i,j,kLev,bi,bj)           gU(i,j,kLev,bi,bj)=gU(i,j,kLev,bi,bj)
65       &                      -kV*uVel(i,j,kLev,bi,bj)       &                     -kV*uVel(i,j,kLev,bi,bj)
66          ENDIF          ENDIF
67         ENDDO         ENDDO
68        ENDDO        ENDDO
69    
70        RETURN        RETURN
71        END        END
72  CStartOfInterface  
73    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
74    CBOP
75    C     !ROUTINE: EXTERNAL_FORCING_V
76    C     !INTERFACE:
77        SUBROUTINE EXTERNAL_FORCING_V(        SUBROUTINE EXTERNAL_FORCING_V(
78       I           iMin, iMax, jMin, jMax,bi,bj,kLev,       I           iMin,iMax, jMin,jMax, bi,bj, kLev,
79       I           myCurrentTime,myThid)       I           myTime, myThid )
80  C     /==========================================================\  C     !DESCRIPTION: \bv
81  C     | S/R EXTERNAL_FORCING_V                                   |  C     *==========================================================*
82  C     | o Contains problem specific forcing for merid velocity.  |  C     | S/R EXTERNAL_FORCING_V
83  C     |==========================================================|  C     | o Contains problem specific forcing for merid velocity.
84  C     | Adds terms to gV for forcing by external sources         |  C     *==========================================================*
85  C     | e.g. wind stress, bottom friction etc..................  |  C     | Adds terms to gV for forcing by external sources
86  C     \==========================================================/  C     | e.g. wind stress, bottom friction etc ...
87        IMPLICIT NONE  C     *==========================================================*
88    C     \ev
89    
90    C     !USES:
91          IMPLICIT NONE
92  C     == Global data ==  C     == Global data ==
93  #include "SIZE.h"  #include "SIZE.h"
94  #include "EEPARAMS.h"  #include "EEPARAMS.h"
# Line 85  C     == Global data == Line 97  C     == Global data ==
97  #include "DYNVARS.h"  #include "DYNVARS.h"
98  #include "FFIELDS.h"  #include "FFIELDS.h"
99    
100    C     !INPUT/OUTPUT PARAMETERS:
101  C     == Routine arguments ==  C     == Routine arguments ==
102  C     iMin - Working range of tile for applying forcing.  C     iMin,iMax :: Working range of x-index for applying forcing.
103  C     iMax  C     jMin,jMax :: Working range of y-index for applying forcing.
104  C     jMin  C     bi,bj     :: Current tile indices
105  C     jMax  C     kLev      :: Current vertical level index
106  C     kLev  C     myTime    :: Current time in simulation
107    C     myThid    :: Thread Id number
108        INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj        INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
109        _RL myCurrentTime        _RL myTime
110        INTEGER myThid        INTEGER myThid
111  CEndOfInterface  
112    C     !LOCAL VARIABLES:
113  C     == Local variables ==  C     == Local variables ==
114  C     Loop counters  C     i,j       :: Loop counters
115        INTEGER I, J        INTEGER i, j
116  C     _RL uKf  CEOP
117  C     _RL levelOfGround        _RL recip_P0g, termP, kV, kF, sigma_b
118  C     _RL criticalLevel  
119  C     _RL levelOfVelPoint  C--   Forcing term(s)
120  C     _RL dist1        kF=1. _d 0/86400. _d 0
121  C     _RL dist2        sigma_b = 0.7 _d 0
122  C     _RL decayFac        DO j=1,sNy+1
123  C     _RL velDragHeightFac  c      DO i=1,sNx
124        _RL termP,kV,kF  C-jmc: Without CD-scheme, this is OK ; but with CD-scheme, needs to cover [0:sNx+1]
125           DO i=0,sNx+1
126        kF=1./86400.          IF ( hFacS(i,j,kLev,bi,bj) .GT. 0. ) THEN
127        DO J=jMin,jMax           recip_P0g=MAX(recip_Rcol(i,j,bi,bj),recip_Rcol(i,j-1,bi,bj))
128         DO I=iMin,iMax           termP=0.5 _d 0*( MIN(rF(kLev)*recip_P0g,1. _d 0)
129          IF ( HFacS(i,j,kLev,bi,bj) .GT. 0. ) THEN       &                   +rF(kLev+1)*recip_P0g )
130  C        termP=0.5*( rF(kLev) + min( rF(kLev+1) ,  c        termP=0.5 _d 0*( rF(kLev) + rF(kLev+1) )*recip_P0g
131  C    &           min(H(I,J,bi,bj),H(I,J-1,bi,bj))            ) )           kV=kF*MAX( 0. _d 0, (termP-sigma_b)/(1. _d 0-sigma_b) )
          termP=0.5*( rF(kLev) + rF(kLev+1) )  
 C        termP=rC(kLev)  
          kV=kF*MAX(0., (termP*recip_H(I,J,bi,bj)-0.7)/(1.-0.7) )  
132           gV(i,j,kLev,bi,bj)=gV(i,j,kLev,bi,bj)           gV(i,j,kLev,bi,bj)=gV(i,j,kLev,bi,bj)
133       &                      -kV*vVel(i,j,kLev,bi,bj)       &                      -kV*vVel(i,j,kLev,bi,bj)
134          ENDIF          ENDIF
# Line 126  C        termP=rC(kLev) Line 137  C        termP=rC(kLev)
137    
138        RETURN        RETURN
139        END        END
140  CStartOfInterface  
141    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
142    CBOP
143    C     !ROUTINE: EXTERNAL_FORCING_T
144    C     !INTERFACE:
145        SUBROUTINE EXTERNAL_FORCING_T(        SUBROUTINE EXTERNAL_FORCING_T(
146       I           iMin, iMax, jMin, jMax,bi,bj,kLev,       I           iMin,iMax, jMin,jMax, bi,bj, kLev,
147       I           maskC,       I           myTime, myThid )
148       I           myCurrentTime,myThid)  C     !DESCRIPTION: \bv
149  C     /==========================================================\  C     *==========================================================*
150  C     | S/R EXTERNAL_FORCING_T                                   |  C     | S/R EXTERNAL_FORCING_T
151  C     | o Contains problem specific forcing for temperature.     |  C     | o Contains problem specific forcing for temperature.
152  C     |==========================================================|  C     *==========================================================*
153  C     | Adds terms to gT for forcing by external sources         |  C     | Adds terms to gT for forcing by external sources
154  C     | e.g. heat flux, climatalogical relaxation..............  |  C     | e.g. heat flux, climatalogical relaxation, etc ...
155  C     \==========================================================/  C     *==========================================================*
156        IMPLICIT NONE  C     \ev
157    
158    C     !USES:
159          IMPLICIT NONE
160  C     == Global data ==  C     == Global data ==
161  #include "SIZE.h"  #include "SIZE.h"
162  #include "EEPARAMS.h"  #include "EEPARAMS.h"
# Line 148  C     == Global data == Line 165  C     == Global data ==
165  #include "DYNVARS.h"  #include "DYNVARS.h"
166  #include "FFIELDS.h"  #include "FFIELDS.h"
167    
168    C     !INPUT/OUTPUT PARAMETERS:
169  C     == Routine arguments ==  C     == Routine arguments ==
170  C     iMin - Working range of tile for applying forcing.  C     iMin,iMax :: Working range of x-index for applying forcing.
171  C     iMax  C     jMin,jMax :: Working range of y-index for applying forcing.
172  C     jMin  C     bi,bj     :: Current tile indices
173  C     jMax  C     kLev      :: Current vertical level index
174  C     kLev  C     myTime    :: Current time in simulation
175        _RS maskC (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  C     myThid    :: Thread Id number
176        INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj        INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
177        _RL myCurrentTime        _RL myTime
178        INTEGER myThid        INTEGER myThid
 CEndOfInterface  
179    
180    C     !LOCAL VARIABLES:
181  C     == Local variables ==  C     == Local variables ==
182  C     Loop counters  C     i,j       :: Loop counters
183        INTEGER I, J        INTEGER i, j
184        _RL thetaLim,kT,ka,ks,term1,term2,thetaEq,termP,rSurf  CEOP
185          _RL thetaLim,kT,ka,ks,sigma_b,term1,term2,thetaEq,termP
186        rSurf=1.E5  
187        ka=1./(40.*86400.)  C--   Forcing term(s)
188        ks=1./(4. *86400.)        ka=1. _d 0/(40. _d 0*86400. _d 0)
189        DO J=jMin,jMax        ks=1. _d 0/(4. _d 0 *86400. _d 0)
190         term1=60.*(sin(yC(1,J,bi,bj)*deg2rad)**2)        sigma_b = 0.7 _d 0
191  C      termP=0.5*( rF(kLev) + min( rF(kLev+1) , H(I,J,bi,bj) ) )        DO j=1,sNy
192         termP=0.5*( rF(kLev) + rF(kLev+1) )         DO i=1,sNx
193  C      termP=rC(kLev)           term1=60. _d 0*(SIN(yC(i,j,bi,bj)*deg2rad)**2)
194         term2=10.*log(termP/rSurf)           termP=0.5 _d 0*( rF(kLev) + rF(kLev+1) )
195       &          *(cos(yC(1,J,bi,bj)*deg2rad)**2)           term2=10. _d 0*LOG(termP/atm_po)
196         thetaLim = 200. / ((termP/rSurf)**(2./7.))       &            *(COS(yC(i,j,bi,bj)*deg2rad)**2)
197         thetaEq=315.-term1-term2           thetaLim = 200. _d 0/ ((termP/atm_po)**atm_kappa)
198         thetaEq=MAX(thetaLim,thetaEq)           thetaEq=315. _d 0-term1-term2
199         DO I=iMin,iMax           thetaEq=MAX(thetaLim,thetaEq)
200          kT=ka+(ks-ka)           termP=0.5 _d 0*( MIN(rF(kLev),Ro_surf(i,j,bi,bj))+rF(kLev+1) )
201       &    *MAX(0., (termP*recip_H(I,J,bi,bj)-0.7)/(1.-0.7) )           kT=ka+(ks-ka)
202       &    *COS((yC(1,J,bi,bj)*deg2rad))**4       &     *MAX(0. _d 0,
203         &       (termP*recip_Rcol(i,j,bi,bj)-sigma_b)/(1. _d 0-sigma_b) )
204         &     *COS((yC(i,j,bi,bj)*deg2rad))**4
205           gT(i,j,kLev,bi,bj)=gT(i,j,kLev,bi,bj)           gT(i,j,kLev,bi,bj)=gT(i,j,kLev,bi,bj)
206       &        - kT*( theta(I,J,kLev,bi,bj)-thetaEq )       &        - kT*( theta(i,j,kLev,bi,bj)-thetaEq )
207       &            *maskC(i,j)       &            *maskC(i,j,kLev,bi,bj)
208         ENDDO         ENDDO
209        ENDDO        ENDDO
210    
211        RETURN        RETURN
212        END        END
213  CStartOfInterface  
214    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
215    CBOP
216    C     !ROUTINE: EXTERNAL_FORCING_S
217    C     !INTERFACE:
218        SUBROUTINE EXTERNAL_FORCING_S(        SUBROUTINE EXTERNAL_FORCING_S(
219       I           iMin, iMax, jMin, jMax,bi,bj,kLev,       I           iMin,iMax, jMin,jMax, bi,bj, kLev,
220       I           maskC,       I           myTime, myThid )
221       I           myCurrentTime,myThid)  
222  C     /==========================================================\  C     !DESCRIPTION: \bv
223  C     | S/R EXTERNAL_FORCING_S                                   |  C     *==========================================================*
224  C     | o Contains problem specific forcing for merid velocity.  |  C     | S/R EXTERNAL_FORCING_S
225  C     |==========================================================|  C     | o Contains problem specific forcing for merid velocity.
226  C     | Adds terms to gS for forcing by external sources         |  C     *==========================================================*
227  C     | e.g. fresh-water flux, climatalogical relaxation.......  |  C     | Adds terms to gS for forcing by external sources
228  C     \==========================================================/  C     | e.g. fresh-water flux, climatalogical relaxation, etc ...
229        IMPLICIT NONE  C     *==========================================================*
230    C     \ev
231    
232    C     !USES:
233          IMPLICIT NONE
234  C     == Global data ==  C     == Global data ==
235  #include "SIZE.h"  #include "SIZE.h"
236  #include "EEPARAMS.h"  #include "EEPARAMS.h"
# Line 212  C     == Global data == Line 239  C     == Global data ==
239  #include "DYNVARS.h"  #include "DYNVARS.h"
240  #include "FFIELDS.h"  #include "FFIELDS.h"
241    
242    C     !INPUT/OUTPUT PARAMETERS:
243  C     == Routine arguments ==  C     == Routine arguments ==
244  C     iMin - Working range of tile for applying forcing.  C     iMin,iMax :: Working range of x-index for applying forcing.
245  C     iMax  C     jMin,jMax :: Working range of y-index for applying forcing.
246  C     jMin  C     bi,bj     :: Current tile indices
247  C     jMax  C     kLev      :: Current vertical level index
248  C     kLev  C     myTime    :: Current time in simulation
249        _RS maskC (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  C     myThid    :: Thread Id number
250        INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj        INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
251        _RL myCurrentTime        _RL myTime
252        INTEGER myThid        INTEGER myThid
 CEndOfInterface  
253    
254    C     !LOCAL VARIABLES:
255  C     == Local variables ==  C     == Local variables ==
256  C     Loop counters  C     i,j       :: Loop counters
257        INTEGER I, J  c     INTEGER i, j
258    CEOP
259    
260    C--   Forcing term(s)
261    
262        RETURN        RETURN
263        END        END
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.8
  ViewVC Help
Powered by ViewVC 1.1.22