/[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.5 by jmc, Fri Jul 6 22:13:37 2001 UTC revision 1.6 by jmc, Sun Jul 17 18:40:10 2005 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"
# Line 24  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
 C     _RL criticalLevel  
 C     _RL levelOfVelPoint  
 C     _RL dist1  
 C     _RL dist2  
 C     _RL decayFac  
 C     _RL velDragHeightFac  
       _RL recip_P0g,termP,kV,kF,sigma_b  
50    
51  C--   Forcing term(s)  C--   Forcing term(s)
52        kF=1. _d 0/86400. _d 0        kF=1. _d 0/86400. _d 0
53        sigma_b = 0.7 _d 0        sigma_b = 0.7 _d 0
54  c     DO J=jMin,jMax  c     DO j=1,sNy
55  c      DO I=iMin,iMax  C-jmc: Without CD-scheme, this is OK ; but with CD-scheme, needs to cover [0:sNy+1]
56        DO J=1,sNy        DO j=0,sNy+1
57         DO I=1,sNx+1         DO i=1,sNx+1
58          IF ( hFacW(i,j,kLev,bi,bj) .GT. 0. ) THEN          IF ( hFacW(i,j,kLev,bi,bj) .GT. 0. ) THEN
59           recip_P0g=MAX(recip_Rcol(I,J,bi,bj),recip_Rcol(I-1,J,bi,bj))           recip_P0g=MAX(recip_Rcol(i,j,bi,bj),recip_Rcol(i-1,j,bi,bj))
60           termP=0.5 _d 0*( MIN(rF(kLev)*recip_P0g,1. _d 0)           termP=0.5 _d 0*( MIN(rF(kLev)*recip_P0g,1. _d 0)
61       &                   +rF(kLev+1)*recip_P0g )       &                   +rF(kLev+1)*recip_P0g )
62  c        termP=0.5 _d 0*( rF(kLev) + rF(kLev+1) )*recip_P0g  c        termP=0.5 _d 0*( rF(kLev) + rF(kLev+1) )*recip_P0g
63           kV=kF*MAX( 0. _d 0, (termP-sigma_b)/(1. _d 0-sigma_b) )           kV=kF*MAX( 0. _d 0, (termP-sigma_b)/(1. _d 0-sigma_b) )
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 90  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
 C     _RL criticalLevel  
 C     _RL levelOfVelPoint  
 C     _RL dist1  
 C     _RL dist2  
 C     _RL decayFac  
 C     _RL velDragHeightFac  
       _RL recip_P0g,termP,kV,kF,sigma_b  
118    
119  C--   Forcing term(s)  C--   Forcing term(s)
120        kF=1. _d 0/86400. _d 0        kF=1. _d 0/86400. _d 0
121        sigma_b = 0.7 _d 0        sigma_b = 0.7 _d 0
122  c     DO J=jMin,jMax        DO j=1,sNy+1
123  c      DO I=iMin,iMax  c      DO i=1,sNx
124        DO J=1,sNy+1  C-jmc: Without CD-scheme, this is OK ; but with CD-scheme, needs to cover [0:sNx+1]
125         DO I=1,sNx         DO i=0,sNx+1
126          IF ( hFacS(i,j,kLev,bi,bj) .GT. 0. ) THEN          IF ( hFacS(i,j,kLev,bi,bj) .GT. 0. ) THEN
127           recip_P0g=MAX(recip_Rcol(I,J,bi,bj),recip_Rcol(I,J-1,bi,bj))           recip_P0g=MAX(recip_Rcol(i,j,bi,bj),recip_Rcol(i,j-1,bi,bj))
128           termP=0.5 _d 0*( MIN(rF(kLev)*recip_P0g,1. _d 0)           termP=0.5 _d 0*( MIN(rF(kLev)*recip_P0g,1. _d 0)
129       &                   +rF(kLev+1)*recip_P0g )       &                   +rF(kLev+1)*recip_P0g )
130  c        termP=0.5 _d 0*( rF(kLev) + rF(kLev+1) )*recip_P0g  c        termP=0.5 _d 0*( rF(kLev) + rF(kLev+1) )*recip_P0g
# Line 135  c        termP=0.5 _d 0*( rF(kLev) + rF( Line 137  c        termP=0.5 _d 0*( rF(kLev) + rF(
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           myCurrentTime,myThid)       I           myTime, myThid )
148  C     /==========================================================\  C     !DESCRIPTION: \bv
149  C     | S/R EXTERNAL_FORCING_T                                   |  C     *==========================================================*
150  C     | o Contains problem specific forcing for temperature.     |  C     | S/R EXTERNAL_FORCING_T
151  C     |==========================================================|  C     | o Contains problem specific forcing for temperature.
152  C     | Adds terms to gT for forcing by external sources         |  C     *==========================================================*
153  C     | e.g. heat flux, climatalogical relaxation..............  |  C     | Adds terms to gT for forcing by external sources
154  C     \==========================================================/  C     | e.g. heat flux, climatalogical relaxation, etc ...
155        IMPLICIT NONE  C     *==========================================================*
156    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 156  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    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    CEOP
185        _RL thetaLim,kT,ka,ks,sigma_b,term1,term2,thetaEq,termP        _RL thetaLim,kT,ka,ks,sigma_b,term1,term2,thetaEq,termP
186    
187  C--   Forcing term(s)  C--   Forcing term(s)
188        ka=1. _d 0/(40. _d 0*86400. _d 0)        ka=1. _d 0/(40. _d 0*86400. _d 0)
189        ks=1. _d 0/(4. _d 0 *86400. _d 0)        ks=1. _d 0/(4. _d 0 *86400. _d 0)
190        sigma_b = 0.7 _d 0        sigma_b = 0.7 _d 0
191        DO J=jMin,jMax        DO j=1,sNy
192         DO I=iMin,iMax         DO i=1,sNx
193           term1=60. _d 0*(sin(yC(I,J,bi,bj)*deg2rad)**2)           term1=60. _d 0*(SIN(yC(i,j,bi,bj)*deg2rad)**2)
194           termP=0.5 _d 0*( rF(kLev) + rF(kLev+1) )           termP=0.5 _d 0*( rF(kLev) + rF(kLev+1) )
195           term2=10. _d 0*log(termP/atm_po)           term2=10. _d 0*LOG(termP/atm_po)
196       &            *(cos(yC(I,J,bi,bj)*deg2rad)**2)       &            *(COS(yC(i,j,bi,bj)*deg2rad)**2)
197           thetaLim = 200. _d 0/ ((termP/atm_po)**atm_kappa)           thetaLim = 200. _d 0/ ((termP/atm_po)**atm_kappa)
198           thetaEq=315. _d 0-term1-term2           thetaEq=315. _d 0-term1-term2
199           thetaEq=MAX(thetaLim,thetaEq)           thetaEq=MAX(thetaLim,thetaEq)
200           termP=0.5 _d 0*( MIN(rF(kLev),Ro_surf(I,J,bi,bj))+rF(kLev+1) )           termP=0.5 _d 0*( MIN(rF(kLev),Ro_surf(i,j,bi,bj))+rF(kLev+1) )
201           kT=ka+(ks-ka)           kT=ka+(ks-ka)
202       &     *MAX(0. _d 0,       &     *MAX(0. _d 0,
203       &       (termP*recip_Rcol(I,J,bi,bj)-sigma_b)/(1. _d 0-sigma_b) )       &       (termP*recip_Rcol(i,j,bi,bj)-sigma_b)/(1. _d 0-sigma_b) )
204       &     *COS((yC(I,J,bi,bj)*deg2rad))**4       &     *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,kLev,bi,bj)       &            *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           myCurrentTime,myThid)       I           myTime, myThid )
221  C     /==========================================================\  
222  C     | S/R EXTERNAL_FORCING_S                                   |  C     !DESCRIPTION: \bv
223  C     | o Contains problem specific forcing for merid velocity.  |  C     *==========================================================*
224  C     |==========================================================|  C     | S/R EXTERNAL_FORCING_S
225  C     | Adds terms to gS for forcing by external sources         |  C     | o Contains problem specific forcing for merid velocity.
226  C     | e.g. fresh-water flux, climatalogical relaxation.......  |  C     *==========================================================*
227  C     \==========================================================/  C     | Adds terms to gS for forcing by external sources
228        IMPLICIT NONE  C     | e.g. fresh-water flux, climatalogical relaxation, etc ...
229    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 219  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    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)  C--   Forcing term(s)
261    
Legend:
Removed from v.1.5  
changed lines
  Added in v.1.6
  ViewVC Help
Powered by ViewVC 1.1.22