/[MITgcm]/MITgcm_contrib/ksnow/press_release/code/update_cg2d.F
ViewVC logotype

Diff of /MITgcm_contrib/ksnow/press_release/code/update_cg2d.F

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

revision 1.1 by ksnow, Fri Dec 16 15:23:18 2016 UTC revision 1.4 by dgoldberg, Sat Mar 4 11:57:16 2017 UTC
# Line 7  C $Name$ Line 7  C $Name$
7  # include "SHELFICE_OPTIONS.h"  # include "SHELFICE_OPTIONS.h"
8  #endif  #endif
9    
   
10  CBOP  CBOP
11  C     !ROUTINE: UPDATE_CG2D  C     !ROUTINE: UPDATE_CG2D
12  C     !INTERFACE:  C     !INTERFACE:
# Line 33  C     === Global variables === Line 32  C     === Global variables ===
32  #include "GRID.h"  #include "GRID.h"
33  #include "SURFACE.h"  #include "SURFACE.h"
34  #include "CG2D.h"  #include "CG2D.h"
35    #ifdef ALLOW_SOLVE4_PS_AND_DRAG
36    # include "DYNVARS.h"
37    #endif /* ALLOW_SOLVE4_PS_AND_DRAG */
38  #ifdef ALLOW_PRESSURE_RELEASE_CODE  #ifdef ALLOW_PRESSURE_RELEASE_CODE
39  # include "SHELFICE.h"  # include "SHELFICE.h"
40  #endif  #endif
 #ifdef IMPLICIT_BOTTOMSIDEDRAG  
 # include "MOM_VISC.h"  
 #endif  
   
41    
42  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
43  C     === Routine arguments ===  C     === Routine arguments ===
# Line 61  C     faceArea :: Temporary used to hold Line 59  C     faceArea :: Temporary used to hold
59        INTEGER i, j, k, ks        INTEGER i, j, k, ks
60        _RL     faceArea        _RL     faceArea
61        _RL     pW_tmp, pS_tmp        _RL     pW_tmp, pS_tmp
62        _RL depth, eff_depth, depth_fac, pinit        LOGICAL updatePreCond
63  #ifdef ALLOW_PRESSURE_RELEASE_CODE  #ifdef ALLOW_PRESSURE_RELEASE_CODE
64        _RL     curPrelVisc, avgGrd        _RL     curPrelVisc, avgGrd
65          _RL depth, eff_depth, depth_fac, pinit
66          _RL drag_fac
67  #endif  #endif
       LOGICAL updatePreCond  
68  CEOP  CEOP
69    
70  C--   Decide when to update cg2d Preconditioner :  C--   Decide when to update cg2d Preconditioner :
# Line 90  C     aS2d: integral in Z Ay/dY Line 89  C     aS2d: integral in Z Ay/dY
89  #endif /* ALLOW_AUTODIFF */  #endif /* ALLOW_AUTODIFF */
90           ENDDO           ENDDO
91          ENDDO          ENDDO
92          DO k=1,Nr  #ifdef ALLOW_SOLVE4_PS_AND_DRAG
93           DO j=1,sNy+1          IF ( selectImplicitDrag.EQ.2 ) THEN
94            DO i=1,sNx+1           DO k=1,Nr
95              DO j=1,sNy+1
96               DO i=1,sNx+1
97                faceArea = _dyG(i,j,bi,bj)*drF(k)*deepFacC(k)*rhoFacC(k)
98         &                *_hFacW(i,j,k,bi,bj)
99    #ifdef ALLOW_PRESSURE_RELEASE_CODE
100                IF (depthColW(i,j,bi,bj).lt.cg2dminColumnEps) THEN        
101    
102                
103                 drag_fac = _recip_hFacW(i,j,k,bi,bj)**2*recip_drF(k)**2*
104         &        pReleaseDamp * viscArNr(k) *
105         &        1./(1+exp(-10./cg2dminColumnEps*
106         &         (-1.)*(depthColW(i,j,bi,bj)-cg2dminColumnEps/2)))
107    
108                 aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj)
109         &               + faceArea*dU_psFacX(i,j,k,bi,bj)
110         &               / (1 + drag_fac*deltaTmom)
111         &                         *recip_dxC(i,j,bi,bj)
112                ELSE
113    #endif
114                 aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj)
115         &               + faceArea*dU_psFacX(i,j,k,bi,bj)
116         &                         *recip_dxC(i,j,bi,bj)
117    #ifdef ALLOW_PRESSURE_RELEASE_CODE
118                ENDIF
119    #endif
120    
121                faceArea = _dxG(i,j,bi,bj)*drF(k)*deepFacC(k)*rhoFacC(k)
122         &                *_hFacS(i,j,k,bi,bj)
123    #ifdef ALLOW_PRESSURE_RELEASE_CODE
124                IF (depthColS(i,j,bi,bj).lt.cg2dminColumnEps) THEN        
125    
126                 drag_fac = _recip_hFacS(i,j,k,bi,bj)**2*recip_drF(k)**2*
127         &        pReleaseDamp * viscArNr(k)*
128         &        1./(1+exp(-10./cg2dminColumnEps*
129         &         (-1.)*(depthColS(i,j,bi,bj)-cg2dminColumnEps/2)))
130    
131                
132                 aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)
133         &               + faceArea*dV_psFacY(i,j,k,bi,bj)
134         &               / (1 + drag_fac*deltaTmom)
135         &                         *recip_dyC(i,j,bi,bj)
136                ELSE
137    #endif
138                 aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)
139         &               + faceArea*dV_psFacY(i,j,k,bi,bj)
140         &                         *recip_dyC(i,j,bi,bj)
141    #ifdef ALLOW_PRESSURE_RELEASE_CODE
142                ENDIF
143    #endif
144               ENDDO
145              ENDDO
146             ENDDO
147            ELSE
148    #endif /* ALLOW_SOLVE4_PS_AND_DRAG */
149             DO k=1,Nr
150              DO j=1,sNy+1
151               DO i=1,sNx+1
152  C  deep-model: *deepFacC (faceArea), /deepFacC (recip_dx,y): => no net effect  C  deep-model: *deepFacC (faceArea), /deepFacC (recip_dx,y): => no net effect
153             faceArea = _dyG(i,j,bi,bj)*drF(k)  
154       &               *_hFacW(i,j,k,bi,bj)  #ifdef ALLOW_PRESSURE_RELEASE_CODE
155             aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj)              IF (depthColW(i,j,bi,bj).lt.cg2dminColumnEps) THEN        
156  #ifndef IMPLICIT_BOTTOMSIDEDRAG               drag_fac = _recip_hFacW(i,j,k,bi,bj)**2*recip_drF(k)**2*
157       &              + faceArea*recip_dxC(i,j,bi,bj)       &        pReleaseDamp * viscArNr(k)*
158  #else       &        1./(1+exp(-10./cg2dminColumnEps*
159       &              + faceArea/       &         (-1.)*(depthColW(i,j,bi,bj)-cg2dminColumnEps/2)))
160       &                (1 - uDragTermsCommon(i,j,k,bi,bj)*deltaTmom) *              ELSE
161       &                recip_dxC(i,j,bi,bj)               drag_fac = 0. _d 0
162  #endif              ENDIF
   
            faceArea = _dxG(i,j,bi,bj)*drF(k)  
      &               *_hFacS(i,j,k,bi,bj)  
            aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)  
 #ifndef IMPLICIT_BOTTOMSIDEDRAG  
      &              + faceArea*recip_dyC(i,j,bi,bj)  
 #else  
      &              + faceArea/  
      &                (1 - vDragTermsCommon(i,j,k,bi,bj)*deltaTmom) *  
      &                recip_dyC(i,j,bi,bj)  
163  #endif  #endif
164                faceArea = _dyG(i,j,bi,bj)*drF(k)
165         &                *_hFacW(i,j,k,bi,bj)
166    
167                
168                aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj)
169         &               + faceArea*recip_dxC(i,j,bi,bj)
170    #ifdef ALLOW_PRESSURE_RELEASE_CODE
171         &       / (1 + drag_fac*deltaTmom)        
172    #endif
173    #ifdef ALLOW_PRESSURE_RELEASE_CODE
174                IF (depthColS(i,j,bi,bj).lt.cg2dminColumnEps) THEN        
175                 drag_fac = _recip_hFacS(i,j,k,bi,bj)**2*recip_drF(k)**2*
176         &        pReleaseDamp * viscArNr(k) *
177         &        1./(1+exp(-10./cg2dminColumnEps*
178         &         (-1.)*(depthColS(i,j,bi,bj)-cg2dminColumnEps/2)))
179                ELSE
180                 drag_fac = 0. _d 0
181                ENDIF
182    #endif
183                faceArea = _dxG(i,j,bi,bj)*drF(k)
184         &                *_hFacS(i,j,k,bi,bj)
185    
186                aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)
187         &               + faceArea*recip_dyC(i,j,bi,bj)
188    #ifdef ALLOW_PRESSURE_RELEASE_CODE
189         &       / (1 + drag_fac*deltaTmom)        
190    #endif
191               ENDDO
192            ENDDO            ENDDO
193           ENDDO           ENDDO
194          ENDDO  #ifdef ALLOW_SOLVE4_PS_AND_DRAG
195            ENDIF
196    #endif /* ALLOW_SOLVE4_PS_AND_DRAG */
197          DO j=1,sNy+1          DO j=1,sNy+1
198           DO i=1,sNx+1           DO i=1,sNx+1
199            aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj)*cg2dNorm            aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj)*cg2dNorm
# Line 125  C  deep-model: *deepFacC (faceArea), /de Line 201  C  deep-model: *deepFacC (faceArea), /de
201  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
202       &                  *maskInC(i,j,bi,bj)*maskInC(i-1,j,bi,bj)       &                  *maskInC(i,j,bi,bj)*maskInC(i-1,j,bi,bj)
203  #endif  #endif
 #ifdef ALLOW_PRESSURE_RELEASE_CODE  
      &                    
 #endif  
204            aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)*cg2dNorm            aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)*cg2dNorm
205       &                     *implicSurfPress*implicDiv2DFlow       &                     *implicSurfPress*implicDiv2DFlow
206  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
# Line 135  C  deep-model: *deepFacC (faceArea), /de Line 208  C  deep-model: *deepFacC (faceArea), /de
208  #endif  #endif
209           ENDDO           ENDDO
210          ENDDO          ENDDO
 !        call write_fld_xy_rl('as2d0','',as2d,0,mythid)  
211  #ifdef ALLOW_PRESSURE_RELEASE_CODE  #ifdef ALLOW_PRESSURE_RELEASE_CODE
212          DO j=1,sNy+1          DO j=1,sNy+1
213           DO i=1,sNx+1           DO i=1,sNx+1
# Line 177  C  deep-model: *deepFacC (faceArea), /de Line 249  C  deep-model: *deepFacC (faceArea), /de
249           ENDDO           ENDDO
250          ENDDO          ENDDO
251  #endif  #endif
 !        call write_fld_xy_rl('as2d0','',as2d,0,mythid)  
252  C--   compute matrix main diagonal :  C--   compute matrix main diagonal :
253          IF ( deepAtmosphere ) THEN          IF ( deepAtmosphere ) THEN
254           DO j=1,sNy           DO j=1,sNy
# Line 187  C--   compute matrix main diagonal : Line 258  C--   compute matrix main diagonal :
258       &       aW2d(i,j,bi,bj) + aW2d(i+1, j ,bi,bj)       &       aW2d(i,j,bi,bj) + aW2d(i+1, j ,bi,bj)
259       &      +aS2d(i,j,bi,bj) + aS2d( i ,j+1,bi,bj)       &      +aS2d(i,j,bi,bj) + aS2d( i ,j+1,bi,bj)
260       &      +freeSurfFac*cg2dNorm*recip_Bo(i,j,bi,bj)*deepFac2F(ks)       &      +freeSurfFac*cg2dNorm*recip_Bo(i,j,bi,bj)*deepFac2F(ks)
261       &                  *rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf       &                  *rA(i,j,bi,bj)/deltaTMom/deltaTFreeSurf
262       &                        )       &                        )
263            ENDDO            ENDDO
264           ENDDO           ENDDO
# Line 198  C--   compute matrix main diagonal : Line 269  C--   compute matrix main diagonal :
269       &       aW2d(i,j,bi,bj) + aW2d(i+1, j ,bi,bj)       &       aW2d(i,j,bi,bj) + aW2d(i+1, j ,bi,bj)
270       &      +aS2d(i,j,bi,bj) + aS2d( i ,j+1,bi,bj)       &      +aS2d(i,j,bi,bj) + aS2d( i ,j+1,bi,bj)
271       &      +freeSurfFac*cg2dNorm*recip_Bo(i,j,bi,bj)       &      +freeSurfFac*cg2dNorm*recip_Bo(i,j,bi,bj)
272       &                  *rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf       &                  *rA(i,j,bi,bj)/deltaTMom/deltaTFreeSurf
273       &                        )       &                        )
274            ENDDO            ENDDO
275           ENDDO           ENDDO
# Line 242  C--   Initialise preconditioner Line 313  C--   Initialise preconditioner
313  C-    if update Preconditioner : end  C-    if update Preconditioner : end
314        ENDIF        ENDIF
315    
   
316        RETURN        RETURN
317        END        END
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.4
  ViewVC Help
Powered by ViewVC 1.1.22