/[MITgcm]/MITgcm/verification/hs94.1x64x5/code/apply_forcing.F
ViewVC logotype

Diff of /MITgcm/verification/hs94.1x64x5/code/apply_forcing.F

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

revision 1.1 by jmc, Fri Jul 11 18:57:32 2014 UTC revision 1.2 by jmc, Wed Aug 20 20:27:43 2014 UTC
# Line 35  C     == Global data == Line 35  C     == Global data ==
35  #include "EEPARAMS.h"  #include "EEPARAMS.h"
36  #include "PARAMS.h"  #include "PARAMS.h"
37  #include "GRID.h"  #include "GRID.h"
38    #include "SURFACE.h"
39  #include "DYNVARS.h"  #include "DYNVARS.h"
40  #include "FFIELDS.h"  #include "FFIELDS.h"
41    
# Line 58  C     !LOCAL VARIABLES: Line 59  C     !LOCAL VARIABLES:
59  C     i,j       :: Loop counters  C     i,j       :: Loop counters
60        INTEGER i, j        INTEGER i, j
61  CEOP  CEOP
62        _RL recip_P0g, termP, kV, kF, sigma_b        _RL recip_P0g, termP, rFullDepth
63          _RL kV, kF, sigma_b
64    
65  C--   Forcing term  C--   Forcing term
66        kF = 1. _d 0/86400. _d 0        kF = 1. _d 0/86400. _d 0
67        sigma_b = 0.7 _d 0        sigma_b = 0.7 _d 0
68          rFullDepth = rF(1)-rF(Nr+1)
69  c     DO j=1,sNy  c     DO j=1,sNy
70  C-jmc: Without CD-scheme, this is OK ; but with CD-scheme, needs to cover [0:sNy+1]  C-jmc: Without CD-scheme, this is OK ; but with CD-scheme, needs to cover [0:sNy+1]
71        DO j=0,sNy+1        DO j=0,sNy+1
72         DO i=1,sNx+1         DO i=1,sNx+1
73          IF ( maskW(i,j,k,bi,bj).EQ.oneRS ) THEN          IF ( maskW(i,j,k,bi,bj).EQ.oneRS ) THEN
74           recip_P0g = MAX(recip_Rcol(i,j,bi,bj),recip_Rcol(i-1,j,bi,bj))           IF ( selectSigmaCoord.EQ.0 ) THEN
75           termP = 0.5 _d 0*( MIN( rF(k)*recip_P0g, oneRL )            recip_P0g = MAX(recip_Rcol(i,j,bi,bj),recip_Rcol(i-1,j,bi,bj))
76       &                     +rF(k+1)*recip_P0g )            termP = 0.5 _d 0*( MIN( rF(k)*recip_P0g, oneRL )
77  c        termP = 0.5 _d 0*( rF(k) + rF(k+1) )*recip_P0g       &                      +rF(k+1)*recip_P0g )
78    c         termP = 0.5 _d 0*( rF(k) + rF(k+1) )*recip_P0g
79             ELSE
80    C-- Pressure at U.point :
81    c         midP = rLowW(i,j,bi,bj) + aHybSigmC(k)*rFullDepth
82    c    &         + bHybSigmC(k)
83    c    &          *(etaHw(i,j,bi,bj)+rSurfW(i,j,bi,bj)-rLowW(i,j,bi,bj))
84    C-- Sigma at U.point :
85    c         termP = ( midP - rLowW(i,j,bi,bj))
86    c    &          /(etaHw(i,j,bi,bj)+rSurfW(i,j,bi,bj)-rLowW(i,j,bi,bj))
87    C-  which simplifies to:
88              termP = aHybSigmC(k)*rFullDepth
89    #ifdef NONLIN_FRSURF
90         &          /(etaHw(i,j,bi,bj)+rSurfW(i,j,bi,bj)-rLowW(i,j,bi,bj))
91    #else
92         &          /(rSurfW(i,j,bi,bj)-rLowW(i,j,bi,bj))
93    #endif
94         &          + bHybSigmC(k)
95             ENDIF
96           kV = kF*MAX( zeroRL, (termP-sigma_b)/(1. _d 0-sigma_b) )           kV = kF*MAX( zeroRL, (termP-sigma_b)/(1. _d 0-sigma_b) )
97           gU_arr(i,j) = gU_arr(i,j)           gU_arr(i,j) = gU_arr(i,j)
98       &               - kV*uVel(i,j,k,bi,bj)       &               - kV*uVel(i,j,k,bi,bj)
# Line 107  C     == Global data == Line 128  C     == Global data ==
128  #include "EEPARAMS.h"  #include "EEPARAMS.h"
129  #include "PARAMS.h"  #include "PARAMS.h"
130  #include "GRID.h"  #include "GRID.h"
131    #include "SURFACE.h"
132  #include "DYNVARS.h"  #include "DYNVARS.h"
133  #include "FFIELDS.h"  #include "FFIELDS.h"
134    
# Line 130  C     !LOCAL VARIABLES: Line 152  C     !LOCAL VARIABLES:
152  C     i,j       :: Loop counters  C     i,j       :: Loop counters
153        INTEGER i, j        INTEGER i, j
154  CEOP  CEOP
155        _RL recip_P0g, termP, kV, kF, sigma_b        _RL recip_P0g, termP, rFullDepth
156          _RL kV, kF, sigma_b
157    
158  C--   Forcing term  C--   Forcing term
159        kF = 1. _d 0/86400. _d 0        kF = 1. _d 0/86400. _d 0
160        sigma_b = 0.7 _d 0        sigma_b = 0.7 _d 0
161          rFullDepth = rF(1)-rF(Nr+1)
162        DO j=1,sNy+1        DO j=1,sNy+1
163  c      DO i=1,sNx  c      DO i=1,sNx
164  C-jmc: Without CD-scheme, this is OK ; but with CD-scheme, needs to cover [0:sNx+1]  C-jmc: Without CD-scheme, this is OK ; but with CD-scheme, needs to cover [0:sNx+1]
165         DO i=0,sNx+1         DO i=0,sNx+1
166          IF ( maskS(i,j,k,bi,bj).EQ.oneRS ) THEN          IF ( maskS(i,j,k,bi,bj).EQ.oneRS ) THEN
167           recip_P0g = MAX(recip_Rcol(i,j,bi,bj),recip_Rcol(i,j-1,bi,bj))           IF ( selectSigmaCoord.EQ.0 ) THEN
168           termP = 0.5 _d 0*( MIN( rF(k)*recip_P0g, oneRL )            recip_P0g = MAX(recip_Rcol(i,j,bi,bj),recip_Rcol(i,j-1,bi,bj))
169       &                     +rF(k+1)*recip_P0g )            termP = 0.5 _d 0*( MIN( rF(k)*recip_P0g, oneRL )
170  c        termP = 0.5 _d 0*( rF(k) + rF(k+1) )*recip_P0g       &                      +rF(k+1)*recip_P0g )
171    c         termP = 0.5 _d 0*( rF(k) + rF(k+1) )*recip_P0g
172             ELSE
173    C-- Pressure at V.point :
174    c         midP = rLowS(i,j,bi,bj) + aHybSigmC(k)*rFullDepth
175    c    &         + bHybSigmC(k)
176    c    &          *(etaHs(i,j,bi,bj)+rSurfS(i,j,bi,bj)-rLowS(i,j,bi,bj))
177    C-- Sigma at V.point :
178    c         termP = ( midP - rLowS(i,j,bi,bj))
179    c    &          /(etaHs(i,j,bi,bj)+rSurfS(i,j,bi,bj)-rLowS(i,j,bi,bj))
180    C-  which simplifies to:
181              termP = aHybSigmC(k)*rFullDepth
182    #ifdef NONLIN_FRSURF
183         &          /(etaHs(i,j,bi,bj)+rSurfS(i,j,bi,bj)-rLowS(i,j,bi,bj))
184    #else
185         &          /(rSurfS(i,j,bi,bj)-rLowS(i,j,bi,bj))
186    #endif
187         &          + bHybSigmC(k)
188             ENDIF
189           kV = kF*MAX( zeroRL, (termP-sigma_b)/(1. _d 0-sigma_b) )           kV = kF*MAX( zeroRL, (termP-sigma_b)/(1. _d 0-sigma_b) )
190           gV_arr(i,j) = gV_arr(i,j)           gV_arr(i,j) = gV_arr(i,j)
191       &               - kV*vVel(i,j,k,bi,bj)       &               - kV*vVel(i,j,k,bi,bj)
# Line 181  C     == Global data == Line 223  C     == Global data ==
223  #include "GRID.h"  #include "GRID.h"
224  #include "DYNVARS.h"  #include "DYNVARS.h"
225  #include "FFIELDS.h"  #include "FFIELDS.h"
 #include "SURFACE.h"  
226    
227  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
228  C     gT_arr    :: the tendency array  C     gT_arr    :: the tendency array
# Line 204  C     i,j       :: Loop counters Line 245  C     i,j       :: Loop counters
245  C     kSurface  :: index of surface level  C     kSurface  :: index of surface level
246        INTEGER i, j        INTEGER i, j
247  CEOP  CEOP
248        _RL thetaLim,kT,ka,ks,sigma_b,term1,term2,thetaEq,termP        _RL thetaLim, kT, ka, ks, sigma_b, term1, term2, thetaEq
249          _RL termP, rFullDepth
250    
251  C--   Forcing term  C--   Forcing term
252        ka = 1. _d 0/(40. _d 0*86400. _d 0)        ka = 1. _d 0/(40. _d 0*86400. _d 0)
253        ks = 1. _d 0/(4. _d 0 *86400. _d 0)        ks = 1. _d 0/(4. _d 0 *86400. _d 0)
254        sigma_b = 0.7 _d 0        sigma_b = 0.7 _d 0
255        DO j=1,sNy        rFullDepth = rF(1)-rF(Nr+1)
256         DO i=1,sNx        DO j=0,sNy+1
257           DO i=0,sNx+1
258           term1 = 60. _d 0*(SIN(yC(i,j,bi,bj)*deg2rad)**2)           term1 = 60. _d 0*(SIN(yC(i,j,bi,bj)*deg2rad)**2)
259           termP = 0.5 _d 0*( rF(k) + rF(k+1) )           termP = 0.5 _d 0*( rF(k) + rF(k+1) )
260           term2 = 10. _d 0*LOG(termP/atm_po)           term2 = 10. _d 0*LOG(termP/atm_po)
261       &            *(COS(yC(i,j,bi,bj)*deg2rad)**2)       &            *(COS(yC(i,j,bi,bj)*deg2rad)**2)
262           thetaLim = 200. _d 0/ ((termP/atm_po)**atm_kappa)           thetaLim = 200. _d 0/ ((termP/atm_po)**atm_kappa)
263           thetaEq = 315. _d 0-term1-term2           thetaEq = 315. _d 0 - term1 - term2
264           thetaEq = MAX(thetaLim,thetaEq)           thetaEq = MAX(thetaLim,thetaEq)
265           termP = 0.5 _d 0*( MIN(rF(k),Ro_surf(i,j,bi,bj))           IF ( selectSigmaCoord.EQ.0 ) THEN
266       &                    + rF(k+1) )            termP = 0.5 _d 0*( MIN(rF(k),Ro_surf(i,j,bi,bj))
267         &                     + rF(k+1) )
268         &                    *recip_Rcol(i,j,bi,bj)
269             ELSE
270    C-- Pressure at T.point :
271    c         midP = R_low(i,j,bi,bj) + aHybSigmC(k)*rFullDepth
272    c    &         + bHybSigmC(k)
273    c    &          *(etaH(i,j,bi,bj)+Ro_surf(i,j,bi,bj)-R_low(i,j,bi,bj))
274    C-- Sigma at T.point :
275    c         termP = ( midP - R_low(i,j,bi,bj))
276    c    &          /(etaH(i,j,bi,bj)+Ro_surf(i,j,bi,bj)-R_low(i,j,bi,bj))
277    C-  which simplifies to:
278              termP = aHybSigmC(k)*rFullDepth
279    #ifdef NONLIN_FRSURF
280         &          /(etaH(i,j,bi,bj)+Ro_surf(i,j,bi,bj)-R_low(i,j,bi,bj))
281    #else
282         &          /(Ro_surf(i,j,bi,bj)-R_low(i,j,bi,bj))
283    #endif
284         &          + bHybSigmC(k)
285             ENDIF
286           kT = ka+(ks-ka)           kT = ka+(ks-ka)
287       &      *MAX( zeroRL,       &      *MAX( zeroRL, (termP-sigma_b)/(1. _d 0-sigma_b) )
      &       (termP*recip_Rcol(i,j,bi,bj)-sigma_b)/(1. _d 0-sigma_b) )  
288       &      *COS((yC(i,j,bi,bj)*deg2rad))**4       &      *COS((yC(i,j,bi,bj)*deg2rad))**4
289           gT_arr(i,j) = gT_arr(i,j)           gT_arr(i,j) = gT_arr(i,j)
290       &               - kT*( theta(i,j,k,bi,bj)-thetaEq )       &               - kT*( theta(i,j,k,bi,bj)-thetaEq )
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
  ViewVC Help
Powered by ViewVC 1.1.22