/[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.2 by ksnow, Mon Jan 30 16:32:17 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         &        30000. * 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                ELSE
112    #endif
113                 aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj)
114         &               + faceArea*dU_psFacX(i,j,k,bi,bj)
115         &                         *recip_dxC(i,j,bi,bj)
116    #ifdef ALLOW_PRESSURE_RELEASE_CODE
117                ENDIF
118    #endif
119    
120                faceArea = _dxG(i,j,bi,bj)*drF(k)*deepFacC(k)*rhoFacC(k)
121         &                *_hFacS(i,j,k,bi,bj)
122    #ifdef ALLOW_PRESSURE_RELEASE_CODE
123                IF (depthColS(i,j,bi,bj).lt.cg2dminColumnEps) THEN        
124    
125                 drag_fac = _recip_hFacS(i,j,k,bi,bj)**2*recip_drF(k)**2*
126         &        30000. * viscArNr(k)*
127         &        1./(1+exp(-10./cg2dminColumnEps*
128         &         (-1.)*(depthColS(i,j,bi,bj)-cg2dminColumnEps/2)))
129    
130                
131                 aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)
132         &               + faceArea*dV_psFacY(i,j,k,bi,bj)
133         &               / (1 + drag_fac*deltaTmom)
134                ELSE
135    #endif
136                 aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)
137         &               + faceArea*dV_psFacY(i,j,k,bi,bj)
138         &                         *recip_dyC(i,j,bi,bj)
139    #ifdef ALLOW_PRESSURE_RELEASE_CODE
140                ENDIF
141    #endif
142               ENDDO
143              ENDDO
144             ENDDO
145            ELSE
146    #endif /* ALLOW_SOLVE4_PS_AND_DRAG */
147             DO k=1,Nr
148              DO j=1,sNy+1
149               DO i=1,sNx+1
150  C  deep-model: *deepFacC (faceArea), /deepFacC (recip_dx,y): => no net effect  C  deep-model: *deepFacC (faceArea), /deepFacC (recip_dx,y): => no net effect
151             faceArea = _dyG(i,j,bi,bj)*drF(k)  
152       &               *_hFacW(i,j,k,bi,bj)  #ifdef ALLOW_PRESSURE_RELEASE_CODE
153             aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj)              IF (depthColW(i,j,bi,bj).lt.cg2dminColumnEps) THEN        
154  #ifndef IMPLICIT_BOTTOMSIDEDRAG               drag_fac = _recip_hFacW(i,j,k,bi,bj)**2*recip_drF(k)**2*
155       &              + faceArea*recip_dxC(i,j,bi,bj)       &        30000. * viscArNr(k)*
156  #else       &        1./(1+exp(-10./cg2dminColumnEps*
157       &              + faceArea/       &         (-1.)*(depthColW(i,j,bi,bj)-cg2dminColumnEps/2)))
158       &                (1 - uDragTermsCommon(i,j,k,bi,bj)*deltaTmom) *              ELSE
159       &                recip_dxC(i,j,bi,bj)               drag_fac = 0. _d 0
160  #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)  
161  #endif  #endif
162                faceArea = _dyG(i,j,bi,bj)*drF(k)
163         &                *_hFacW(i,j,k,bi,bj)
164    
165                
166                aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj)
167         &               + faceArea*recip_dxC(i,j,bi,bj)
168    #ifdef ALLOW_PRESSURE_RELEASE_CODE
169         &       / (1 + drag_fac*deltaTmom)        
170    #endif
171    #ifdef ALLOW_PRESSURE_RELEASE_CODE
172                IF (depthColS(i,j,bi,bj).lt.cg2dminColumnEps) THEN        
173                 drag_fac = _recip_hFacS(i,j,k,bi,bj)**2*recip_drF(k)**2*
174         &        30000. * viscArNr(k) *
175         &        1./(1+exp(-10./cg2dminColumnEps*
176         &         (-1.)*(depthColS(i,j,bi,bj)-cg2dminColumnEps/2)))
177                ELSE
178                 drag_fac = 0. _d 0
179                ENDIF
180    #endif
181                faceArea = _dxG(i,j,bi,bj)*drF(k)
182         &                *_hFacS(i,j,k,bi,bj)
183    
184                aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)
185         &               + faceArea*recip_dyC(i,j,bi,bj)
186    #ifdef ALLOW_PRESSURE_RELEASE_CODE
187         &       / (1 + drag_fac*deltaTmom)        
188    #endif
189               ENDDO
190            ENDDO            ENDDO
191           ENDDO           ENDDO
192          ENDDO  #ifdef ALLOW_SOLVE4_PS_AND_DRAG
193            ENDIF
194    #endif /* ALLOW_SOLVE4_PS_AND_DRAG */
195          DO j=1,sNy+1          DO j=1,sNy+1
196           DO i=1,sNx+1           DO i=1,sNx+1
197            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 199  C  deep-model: *deepFacC (faceArea), /de
199  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
200       &                  *maskInC(i,j,bi,bj)*maskInC(i-1,j,bi,bj)       &                  *maskInC(i,j,bi,bj)*maskInC(i-1,j,bi,bj)
201  #endif  #endif
 #ifdef ALLOW_PRESSURE_RELEASE_CODE  
      &                    
 #endif  
202            aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)*cg2dNorm            aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)*cg2dNorm
203       &                     *implicSurfPress*implicDiv2DFlow       &                     *implicSurfPress*implicDiv2DFlow
204  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
# Line 135  C  deep-model: *deepFacC (faceArea), /de Line 206  C  deep-model: *deepFacC (faceArea), /de
206  #endif  #endif
207           ENDDO           ENDDO
208          ENDDO          ENDDO
 !        call write_fld_xy_rl('as2d0','',as2d,0,mythid)  
209  #ifdef ALLOW_PRESSURE_RELEASE_CODE  #ifdef ALLOW_PRESSURE_RELEASE_CODE
210          DO j=1,sNy+1          DO j=1,sNy+1
211           DO i=1,sNx+1           DO i=1,sNx+1
# Line 177  C  deep-model: *deepFacC (faceArea), /de Line 247  C  deep-model: *deepFacC (faceArea), /de
247           ENDDO           ENDDO
248          ENDDO          ENDDO
249  #endif  #endif
 !        call write_fld_xy_rl('as2d0','',as2d,0,mythid)  
250  C--   compute matrix main diagonal :  C--   compute matrix main diagonal :
251          IF ( deepAtmosphere ) THEN          IF ( deepAtmosphere ) THEN
252           DO j=1,sNy           DO j=1,sNy
# Line 187  C--   compute matrix main diagonal : Line 256  C--   compute matrix main diagonal :
256       &       aW2d(i,j,bi,bj) + aW2d(i+1, j ,bi,bj)       &       aW2d(i,j,bi,bj) + aW2d(i+1, j ,bi,bj)
257       &      +aS2d(i,j,bi,bj) + aS2d( i ,j+1,bi,bj)       &      +aS2d(i,j,bi,bj) + aS2d( i ,j+1,bi,bj)
258       &      +freeSurfFac*cg2dNorm*recip_Bo(i,j,bi,bj)*deepFac2F(ks)       &      +freeSurfFac*cg2dNorm*recip_Bo(i,j,bi,bj)*deepFac2F(ks)
259       &                  *rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf       &                  *rA(i,j,bi,bj)/deltaTMom/deltaTFreeSurf
260       &                        )       &                        )
261            ENDDO            ENDDO
262           ENDDO           ENDDO
# Line 198  C--   compute matrix main diagonal : Line 267  C--   compute matrix main diagonal :
267       &       aW2d(i,j,bi,bj) + aW2d(i+1, j ,bi,bj)       &       aW2d(i,j,bi,bj) + aW2d(i+1, j ,bi,bj)
268       &      +aS2d(i,j,bi,bj) + aS2d( i ,j+1,bi,bj)       &      +aS2d(i,j,bi,bj) + aS2d( i ,j+1,bi,bj)
269       &      +freeSurfFac*cg2dNorm*recip_Bo(i,j,bi,bj)       &      +freeSurfFac*cg2dNorm*recip_Bo(i,j,bi,bj)
270       &                  *rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf       &                  *rA(i,j,bi,bj)/deltaTMom/deltaTFreeSurf
271       &                        )       &                        )
272            ENDDO            ENDDO
273           ENDDO           ENDDO
# Line 242  C--   Initialise preconditioner Line 311  C--   Initialise preconditioner
311  C-    if update Preconditioner : end  C-    if update Preconditioner : end
312        ENDIF        ENDIF
313    
   
314        RETURN        RETURN
315        END        END
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
  ViewVC Help
Powered by ViewVC 1.1.22