/[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.6 by jmc, Sun Jul 17 18:40:10 2005 UTC revision 1.7 by jmc, Fri Sep 24 20:43:35 2010 UTC
# Line 26  C     == Global data == Line 26  C     == Global data ==
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    
# Line 46  C     == Local variables == Line 47  C     == Local variables ==
47  C     i,j       :: Loop counters  C     i,j       :: Loop counters
48        INTEGER i, j        INTEGER i, j
49  CEOP  CEOP
50        _RL recip_P0g, termP, kV, kF, sigma_b        _RL recip_P0g, termP, rFullDepth
51          _RL kV, kF, sigma_b
52    
53  C--   Forcing term(s)  C--   Forcing term(s)
54        kF=1. _d 0/86400. _d 0        kF=1. _d 0/86400. _d 0
55        sigma_b = 0.7 _d 0        sigma_b = 0.7 _d 0
56          rFullDepth = rF(1)-rF(Nr+1)
57  c     DO j=1,sNy  c     DO j=1,sNy
58  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]
59        DO j=0,sNy+1        DO j=0,sNy+1
60         DO i=1,sNx+1         DO i=1,sNx+1
61          IF ( hFacW(i,j,kLev,bi,bj) .GT. 0. ) THEN          IF ( hFacW(i,j,kLev,bi,bj) .GT. 0. ) THEN
62           recip_P0g=MAX(recip_Rcol(i,j,bi,bj),recip_Rcol(i-1,j,bi,bj))           IF ( selectSigmaCoord.EQ.0 ) THEN
63           termP=0.5 _d 0*( MIN(rF(kLev)*recip_P0g,1. _d 0)            recip_P0g=MAX(recip_Rcol(i,j,bi,bj),recip_Rcol(i-1,j,bi,bj))
64       &                   +rF(kLev+1)*recip_P0g )            termP=0.5 _d 0*( MIN(rF(kLev)*recip_P0g,1. _d 0)
65  c        termP=0.5 _d 0*( rF(kLev) + rF(kLev+1) )*recip_P0g       &                    +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) )           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)
# Line 94  C     == Global data == Line 110  C     == Global data ==
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    
# Line 114  C     == Local variables == Line 131  C     == Local variables ==
131  C     i,j       :: Loop counters  C     i,j       :: Loop counters
132        INTEGER i, j        INTEGER i, j
133  CEOP  CEOP
134        _RL recip_P0g, termP, kV, kF, sigma_b        _RL recip_P0g, termP, rFullDepth
135          _RL kV, kF, sigma_b
136    
137  C--   Forcing term(s)  C--   Forcing term(s)
138        kF=1. _d 0/86400. _d 0        kF=1. _d 0/86400. _d 0
139        sigma_b = 0.7 _d 0        sigma_b = 0.7 _d 0
140          rFullDepth = rF(1)-rF(Nr+1)
141        DO j=1,sNy+1        DO j=1,sNy+1
142  c      DO i=1,sNx  c      DO i=1,sNx
143  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]
144         DO i=0,sNx+1         DO i=0,sNx+1
145          IF ( hFacS(i,j,kLev,bi,bj) .GT. 0. ) THEN          IF ( hFacS(i,j,kLev,bi,bj) .GT. 0. ) THEN
146           recip_P0g=MAX(recip_Rcol(i,j,bi,bj),recip_Rcol(i,j-1,bi,bj))           IF ( selectSigmaCoord.EQ.0 ) THEN
147           termP=0.5 _d 0*( MIN(rF(kLev)*recip_P0g,1. _d 0)            recip_P0g=MAX(recip_Rcol(i,j,bi,bj),recip_Rcol(i,j-1,bi,bj))
148       &                   +rF(kLev+1)*recip_P0g )            termP=0.5 _d 0*( MIN(rF(kLev)*recip_P0g,1. _d 0)
149  c        termP=0.5 _d 0*( rF(kLev) + rF(kLev+1) )*recip_P0g       &                    +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) )           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)
# Line 182  C     == Local variables == Line 214  C     == Local variables ==
214  C     i,j       :: Loop counters  C     i,j       :: Loop counters
215        INTEGER i, j        INTEGER i, j
216  CEOP  CEOP
217        _RL thetaLim,kT,ka,ks,sigma_b,term1,term2,thetaEq,termP        _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        ka=1. _d 0/(40. _d 0*86400. _d 0)        ka=1. _d 0/(40. _d 0*86400. _d 0)
222        ks=1. _d 0/(4. _d 0 *86400. _d 0)        ks=1. _d 0/(4. _d 0 *86400. _d 0)
223        sigma_b = 0.7 _d 0        sigma_b = 0.7 _d 0
224          rFullDepth = rF(1)-rF(Nr+1)
225        DO j=1,sNy        DO j=1,sNy
226         DO i=1,sNx         DO i=1,sNx
227           term1=60. _d 0*(SIN(yC(i,j,bi,bj)*deg2rad)**2)           term1=60. _d 0*(SIN(yC(i,j,bi,bj)*deg2rad)**2)
# Line 197  C--   Forcing term(s) Line 231  C--   Forcing term(s)
231           thetaLim = 200. _d 0/ ((termP/atm_po)**atm_kappa)           thetaLim = 200. _d 0/ ((termP/atm_po)**atm_kappa)
232           thetaEq=315. _d 0-term1-term2           thetaEq=315. _d 0-term1-term2
233           thetaEq=MAX(thetaLim,thetaEq)           thetaEq=MAX(thetaLim,thetaEq)
234           termP=0.5 _d 0*( MIN(rF(kLev),Ro_surf(i,j,bi,bj))+rF(kLev+1) )           IF ( selectSigmaCoord.EQ.0 ) THEN
235              termP=0.5 _d 0*( MIN(rF(kLev),Ro_surf(i,j,bi,bj))+rF(kLev+1) )
236         &                  *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)           kT=ka+(ks-ka)
251       &     *MAX(0. _d 0,       &     *MAX(0. _d 0, (termP-sigma_b)/(1. _d 0-sigma_b) )
      &       (termP*recip_Rcol(i,j,bi,bj)-sigma_b)/(1. _d 0-sigma_b) )  
252       &     *COS((yC(i,j,bi,bj)*deg2rad))**4       &     *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 )
Legend:
Removed from v.1.6  
changed lines
  Added in v.1.7
  ViewVC Help
Powered by ViewVC 1.1.22