/[MITgcm]/MITgcm/model/src/ini_masks_etc.F
ViewVC logotype

Diff of /MITgcm/model/src/ini_masks_etc.F

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

revision 1.20 by cnh, Sun Feb 4 14:38:47 2001 UTC revision 1.21 by adcroft, Tue May 29 14:01:37 2001 UTC
# Line 22  C     === Global variables === Line 22  C     === Global variables ===
22  #include "EEPARAMS.h"  #include "EEPARAMS.h"
23  #include "PARAMS.h"  #include "PARAMS.h"
24  #include "GRID.h"  #include "GRID.h"
25    #include "SURFACE.h"
26    
27  C     == Routine arguments ==  C     == Routine arguments ==
28  C     myThid -  Number of this instance of INI_MASKS_ETC  C     myThid -  Number of this instance of INI_MASKS_ETC
# Line 33  C     bi,bj  - Loop counters Line 34  C     bi,bj  - Loop counters
34  C     I,J,K  C     I,J,K
35        INTEGER bi, bj        INTEGER bi, bj
36        INTEGER  I, J, K        INTEGER  I, J, K
       INTEGER kadj_rf, klev_noH  
37  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
38        INTEGER Km1        INTEGER Km1
39        _RL hFacUpper,hFacLower        _RL hFacUpper,hFacLower
40  #endif  #endif
41        _RL hFacCtmp,topo_rkfac        _RL hFacCtmp
42        _RL hFacMnSz        _RL hFacMnSz
43          _RS tmpfld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
44    
45        IF (groundAtK1) THEN  C- Calculate lopping factor hFacC : over-estimate the part inside of the domain
46         topo_rkfac = -rkFac  C    taking into account the lower_R Boundary (Bathymetrie / Top of Atmos)
        kadj_rf = 1  
        klev_noH = Nr+1  
       ELSE  
        topo_rkfac = rkFac  
        kadj_rf = 0  
        klev_noH = 1  
       ENDIF  
   
 C     Calculate lopping factor hFacC  
47        DO bj=myByLo(myThid), myByHi(myThid)        DO bj=myByLo(myThid), myByHi(myThid)
48         DO bi=myBxLo(myThid), myBxHi(myThid)         DO bi=myBxLo(myThid), myBxHi(myThid)
49          DO K=1, Nr          DO K=1, Nr
50             hFacMnSz=max( hFacMin, min(hFacMinDr*recip_drF(k),1. _d 0) )
51           DO J=1-Oly,sNy+Oly           DO J=1-Oly,sNy+Oly
52            DO I=1-Olx,sNx+Olx            DO I=1-Olx,sNx+Olx
53  c          IF (groundAtK1) THEN  C      o Non-dimensional distance between grid bound. and domain lower_R bound.
54  C          o Non-dimensional distance between grid boundary and model depth             hFacCtmp = (rF(K)-R_low(I,J,bi,bj))*recip_drF(K)
55  C            for case with "ground" at K=1 (i.e. like a good atmos model)  C      o Select between, closed, open or partial (0,1,0-1)
 C          e.g. rkfac=+1 => dR/dk<0 (eg. R=p): hFacCtmp=(H(x,y)-rF(k))/drF(K)  
 C          e.g. rkfac=-1 => dR/dk>0 (eg. R=z): hFacCtmp=(rF(K)-H(x,y))/drF(K)  
 c          hFacCtmp=rkFac*(H(I,J,bi,bj)-rF(K+1))*recip_drF(K)  
 c          ELSE  
 C          o Non-dimensional distance between grid boundary and model depth  
 C            for case with "ground" at K=Nr (i.e. like original ocean model)  
 C          e.g. rkfac=+1 => dR/dk<0 (eg. R=z): hFacCtmp=(rF(K)-H(x,y))/drF(K)  
 C          e.g. rkfac=-1 => dR/dk>0 (eg. R=p): hFacCtmp=(H(x,y)-rF(k))/drF(K)  
 c          hFacCtmp=rkFac*(rF(K)-H(I,J,bi,bj))*recip_drF(K)  
 c          ENDIF  
            hFacCtmp=topo_rkfac*(rF(K+kadj_rf)-H(I,J,bi,bj))*recip_drF(K)  
 C          o Select between, closed, open or partial (0,1,0-1)  
56             hFacCtmp=min( max( hFacCtmp, 0. _d 0) , 1. _d 0)             hFacCtmp=min( max( hFacCtmp, 0. _d 0) , 1. _d 0)
57  C          o And there we have it, the fractional open cell volume  C      o Impose minimum fraction and/or size (dimensional)
58             hFacC(I,J,K,bi,bj)=hFacCtmp             IF (hFacCtmp.LT.hFacMnSz) THEN
59  C          o Impose minimum fraction and/or size (dimensional)              IF (hFacCtmp.LT.hFacMnSz*0.5) THEN
            hFacMnSz=max( hFacMin , min(hFacMinDr*recip_drF(k),1. _d 0) )  
            IF (hFacC(I,J,K,bi,bj).LT.hFacMnSz) THEN  
             IF (hFacC(I,J,K,bi,bj).LT.hFacMnSz*0.5) THEN  
60               hFacC(I,J,K,bi,bj)=0.               hFacC(I,J,K,bi,bj)=0.
61              ELSE              ELSE
62               hFacC(I,J,K,bi,bj)=hFacMnSz               hFacC(I,J,K,bi,bj)=hFacMnSz
63              ENDIF              ENDIF
64               ELSE
65                 hFacC(I,J,K,bi,bj)=hFacCtmp
66             ENDIF             ENDIF
            IF (hFacC(I,J,K,bi,bj).NE.0.)  
      &         depthInK(i,j,bi,bj) = depthInK(i,j,bi,bj) + 1.  
67            ENDDO            ENDDO
68           ENDDO           ENDDO
69          ENDDO          ENDDO
70    
71    C-  Re-calculate lower-R Boundary position, taking into account hFacC
72            DO J=1-Oly,sNy+Oly
73             DO I=1-Olx,sNx+Olx
74              R_low(I,J,bi,bj) = rF(1)
75              DO K=Nr,1,-1
76               R_low(I,J,bi,bj) = R_low(I,J,bi,bj)
77         &                      - drF(k)*hFacC(I,J,K,bi,bj)
78              ENDDO
79             ENDDO
80            ENDDO
81    C - end bi,bj loops.
82         ENDDO         ENDDO
83        ENDDO        ENDDO
 C     _EXCH_XYZ_R4(hFacC , myThid )  
 C     _EXCH_XY_R4( depthInK, myThid )  
84    
85        CALL PLOT_FIELD_XYRS( depthInK,  C-  Calculate lopping factor hFacC : Remove part outside of the domain
86       & 'Model Depths K Index' , 1, myThid )  C    taking into account the Reference (=at rest) Surface Position Ro_surf
   
 C Re-calculate depth of ocean, taking into account hFacC  
87        DO bj=myByLo(myThid), myByHi(myThid)        DO bj=myByLo(myThid), myByHi(myThid)
88         DO bi=myBxLo(myThid), myBxHi(myThid)         DO bi=myBxLo(myThid), myBxHi(myThid)
89            DO K=1, Nr
90             hFacMnSz=max( hFacMin, min(hFacMinDr*recip_drF(k),1. _d 0) )
91             DO J=1-Oly,sNy+Oly
92              DO I=1-Olx,sNx+Olx
93    C      o Non-dimensional distance between grid boundary and model surface
94               hFacCtmp = (rF(k)-Ro_surf(I,J,bi,bj))*recip_drF(K)
95    C      o Reduce the previous fraction : substract the outside part.
96               hFacCtmp = hFacC(I,J,K,bi,bj) - max( hFacCtmp, 0. _d 0)
97    C      o set to zero if empty Column :
98               hFacCtmp = max( hFacCtmp, 0. _d 0)
99    C      o Impose minimum fraction and/or size (dimensional)
100               IF (hFacCtmp.LT.hFacMnSz) THEN
101                IF (hFacCtmp.LT.hFacMnSz*0.5) THEN
102                 hFacC(I,J,K,bi,bj)=0.
103                ELSE
104                 hFacC(I,J,K,bi,bj)=hFacMnSz
105                ENDIF
106               ELSE
107                 hFacC(I,J,K,bi,bj)=hFacCtmp
108               ENDIF
109              ENDDO
110             ENDDO
111            ENDDO
112    
113    C-  Re-calculate Reference surface position, taking into account hFacC
114    C   initialize Total column fluid thickness and surface k index
115          DO J=1-Oly,sNy+Oly          DO J=1-Oly,sNy+Oly
116           DO I=1-Olx,sNx+Olx           DO I=1-Olx,sNx+Olx
117            H(I,J,bi,bj)=rF(klev_noH)            tmpfld(I,J,bi,bj) = 0.
118            DO K=1,Nr            k_surf(I,J,bi,bj) = Nr
119             H(I,J,bi,bj)=H(I,J,bi,bj)-            Ro_surf(I,J,bi,bj) = R_low(I,J,bi,bj)
120       &       topo_rkFac*drF(k)*hFacC(I,J,K,bi,bj)            DO K=Nr,1,-1
121               Ro_surf(I,J,bi,bj) = Ro_surf(I,J,bi,bj)
122         &                        + drF(k)*hFacC(I,J,K,bi,bj)
123               IF (hFacC(I,J,K,bi,bj).NE.0.) THEN
124                k_surf(I,J,bi,bj) = k
125                tmpfld(i,j,bi,bj) = tmpfld(i,j,bi,bj) + 1.
126               ENDIF
127            ENDDO            ENDDO
128           ENDDO           ENDDO
129          ENDDO          ENDDO
130    C - end bi,bj loops.
131         ENDDO         ENDDO
132        ENDDO        ENDDO
133  C     _EXCH_XY_R4(H , myThid )  
134        CALL PLOT_FIELD_XYRS(H,'Model depths (ini_masks_etc)',1,myThid)  C     CALL PLOT_FIELD_XYRS( tmpfld,
135        CALL WRITE_FLD_XY_RS( 'Depth',' ',H,0,myThid)  C    &         'Model Depths K Index' , 1, myThid )
136        CALL WRITE_FLD_XYZ_RS( 'hFacC',' ',hFacC,0,myThid)        CALL PLOT_FIELD_XYRS(R_low,
137  C     CALL MDSWRITEFIELD( 'Depth', writeBinaryPrec, .TRUE.,       &         'Model R_low (ini_masks_etc)', 1, myThid)
138  C    &                    'RS', 1, H, 1, -1, myThid )        CALL PLOT_FIELD_XYRS(Ro_surf,
139         &         'Model Ro_surf (ini_masks_etc)', 1, myThid)
140    
141  C     Calculate quantities derived from XY depth map  C     Calculate quantities derived from XY depth map
142        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
143         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
144          DO J=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
145           DO I=1-Olx,sNx+Olx           DO i=1-Olx,sNx+Olx
146  C         Inverse of depth  C         Total fluid column thickness (r_unit) :
147            IF ( h(i,j,bi,bj) .EQ. 0. ) THEN  c         Rcolumn(i,j,bi,bj)= Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
148             recip_H(i,j,bi,bj) = 0.            tmpfld(i,j,bi,bj) = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
149    C         Inverse of fluid column thickness (1/r_unit)
150              IF ( tmpfld(i,j,bi,bj) .LE. 0. ) THEN
151               recip_Rcol(i,j,bi,bj) = 0.
152            ELSE            ELSE
153             recip_H(i,j,bi,bj) = 1. / abs( H(i,j,bi,bj) )             recip_Rcol(i,j,bi,bj) = 1. / tmpfld(i,j,bi,bj)
154            ENDIF            ENDIF
           depthInK(i,j,bi,bj) = 0.  
155           ENDDO           ENDDO
156          ENDDO          ENDDO
157         ENDDO         ENDDO
158        ENDDO        ENDDO
159  C     _EXCH_XY_R4(   recip_H, myThid )  C     _EXCH_XY_R4(   recip_Rcol, myThid )
160          CALL WRITE_FLD_XY_RS( 'Depth',' ',tmpfld,0,myThid)
161          CALL WRITE_FLD_XYZ_RS( 'hFacC',' ',hFacC,0,myThid)
162    C     CALL MDSWRITEFIELD( 'Depth', writeBinaryPrec, .TRUE.,
163    C    &                    'RS', 1, tmpfld, 1, -1, myThid )
164    
165  C     hFacW and hFacS (at U and V points)  C     hFacW and hFacS (at U and V points)
166        DO bj=myByLo(myThid), myByHi(myThid)        DO bj=myByLo(myThid), myByHi(myThid)
# Line 151  C     hFacW and hFacS (at U and V points Line 177  C     hFacW and hFacS (at U and V points
177          ENDDO          ENDDO
178         ENDDO         ENDDO
179        ENDDO        ENDDO
180        _EXCH_XYZ_R4(hFacW , myThid )        CALL EXCH_UV_XYZ_RS(hFacW,hFacS,.FALSE.,myThid)
181        _EXCH_XYZ_R4(hFacS , myThid )  C     The following block allows thin walls representation of non-periodic
182  C     Re-do hFacW and hFacS (at U and V points)  C     boundaries such as happen on the lat-lon grid at the N/S poles.
183    C     We should really supply a flag for doing this.
184        DO bj=myByLo(myThid), myByHi(myThid)        DO bj=myByLo(myThid), myByHi(myThid)
185         DO bi=myBxLo(myThid), myBxHi(myThid)         DO bi=myBxLo(myThid), myBxHi(myThid)
186          DO K=1, Nr          DO K=1, Nr
187           DO J=1-Oly,sNy+Oly           DO J=1-Oly,sNy+Oly
188            DO I=1-Olx+1,sNx+Olx ! Note range            DO I=1-Olx,sNx+Olx
189             hFacW(I,J,K,bi,bj)=             IF (DYG(I,J,bi,bj).EQ.0.) hFacW(I,J,K,bi,bj)=0.
190       &       MIN(hFacC(I,J,K,bi,bj),hFacC(I-1,J,K,bi,bj))             IF (DXG(I,J,bi,bj).EQ.0.) hFacS(I,J,K,bi,bj)=0.
           ENDDO  
          ENDDO  
          DO I=1-Olx,sNx+Olx  
           DO J=1-Oly+1,sNy+Oly ! Note range  
            hFacS(I,J,K,bi,bj)=  
      &       MIN(hFacC(I,J,K,bi,bj),hFacC(I,J-1,K,bi,bj))  
191            ENDDO            ENDDO
192           ENDDO           ENDDO
193          ENDDO          ENDDO
# Line 185  C     Masks and reciprocals of hFac[CWS] Line 206  C     Masks and reciprocals of hFac[CWS]
206            DO I=1-Olx,sNx+Olx            DO I=1-Olx,sNx+Olx
207             IF (HFacC(I,J,K,bi,bj) .NE. 0. ) THEN             IF (HFacC(I,J,K,bi,bj) .NE. 0. ) THEN
208              recip_HFacC(I,J,K,bi,bj) = 1. / HFacC(I,J,K,bi,bj)              recip_HFacC(I,J,K,bi,bj) = 1. / HFacC(I,J,K,bi,bj)
209                maskC(I,J,K,bi,bj) = 1.
210             ELSE             ELSE
211              recip_HFacC(I,J,K,bi,bj) = 0.              recip_HFacC(I,J,K,bi,bj) = 0.
212                maskC(I,J,K,bi,bj) = 0.
213             ENDIF             ENDIF
214             IF (HFacW(I,J,K,bi,bj) .NE. 0. ) THEN             IF (HFacW(I,J,K,bi,bj) .NE. 0. ) THEN
215              recip_HFacW(I,J,K,bi,bj) = 1. / HFacW(I,J,K,bi,bj)              recip_HFacW(I,J,K,bi,bj) = 1. / HFacW(I,J,K,bi,bj)
Legend:
Removed from v.1.20  
changed lines
  Added in v.1.21
  ViewVC Help
Powered by ViewVC 1.1.22