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

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

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

revision 1.4 by jmc, Tue Dec 10 02:55:47 2002 UTC revision 1.5 by jmc, Wed Feb 26 03:11:32 2003 UTC
# Line 31  C     == Global variables == Line 31  C     == Global variables ==
31    
32  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
33  C     == Routine arguments ==  C     == Routine arguments ==
34  C     selectMode :: 0 = find Pfld from Hfld using Theta_Ref profile  C     selectMode ::  > 0 = find Pfld from Hfld ; < 0 = compute Hfld from Pfld
35  C                   1 = find Pfld from Hfld using Analytic fct Theta(Lat,P)  C                   selectFindRoSurf = 0 : use Theta_Ref profile
36  C                  -1 = compute Hfld from Pfld using Analytic  Theta(Lat,P)  C                   selectFindRoSurf = 1 : use analytic fct Theta(Lat,P)
37  C     Hfld (input/outp) :: Ground elevation [m]  C     Hfld (input/outp) :: Ground elevation [m]
38  C     Pfld (outp/input) :: reference Pressure at the ground [Pa]  C     Pfld (outp/input) :: reference Pressure at the ground [Pa]
39        INTEGER selectMode        INTEGER selectMode
# Line 62  C---+----1----+----2----+----3----+----4 Line 62  C---+----1----+----2----+----3----+----4
62        INTEGER bi,bj,i,j,K, Ks        INTEGER bi,bj,i,j,K, Ks
63        _RL ddPI, Po_surf        _RL ddPI, Po_surf
64        _RL phiRef(2*Nr+1), rHalf(2*Nr+1)        _RL phiRef(2*Nr+1), rHalf(2*Nr+1)
65          LOGICAL findPoSurf
66    
67        INTEGER nLevHvR        INTEGER nLevHvR
68        PARAMETER ( nLevHvR = 60 )        PARAMETER ( nLevHvR = 60 )
69        _RL plowHvR, dpHvR, pLevHvR(nLevHvR+1), pMidHvR(nLevHvR)        _RL plowHvR, dpHvR, pLevHvR(nLevHvR+1), pMidHvR(nLevHvR)
70        _RL thetaHvR(nLevHvR), PiHvR(nLevHvR+1), dPiHvR(nLevHvR)        _RL thetaHvR(nLevHvR), PiHvR(nLevHvR+1), dPiHvR(nLevHvR)
71        _RL recip_kappa, PiLoc, zLoc, dzLoc, yLatLoc, phiLoc        _RL recip_kappa, PiLoc, zLoc, dzLoc, yLatLoc, phiLoc
72          _RL  psNorm, rMidKp1
73          _RL ratioRm(Nr), ratioRp(Nr)
74        INTEGER kLev        INTEGER kLev
75  #ifdef CHECK_ANALYLIC_THETA  #ifdef CHECK_ANALYLIC_THETA
76        _RL tmpVar(nLevAn,61)        _RL tmpVar(nLevHvR,61)
77  #endif  #endif
78  CEOP  CEOP
79    
# Line 88  CEOP Line 91  CEOP
91        ENDDO        ENDDO
92         rHalf(2*Nr+1) = rF(Nr+1)         rHalf(2*Nr+1) = rF(Nr+1)
93    
94        IF (selectMode .LE. 0) THEN        IF (selectMode*selectFindRoSurf .LE. 0) THEN
95  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
96  C- Compute Reference Geopotential at Half levels :  C- Compute Reference Geopotential at Half levels :
97  C    Tracer level: phiRef(2K)  ;  Interface_W level: phiRef(2K+1)  C    Tracer level: phiRef(2K)  ;  Interface_W level: phiRef(2K+1)
# Line 143  C- Write to check : Line 146  C- Write to check :
146        ENDDO        ENDDO
147        _END_MASTER( myThid )        _END_MASTER( myThid )
148    
149  C-- endif selectMode =< 0  C-- endif selectMode*selectFindRoSurf =< 0
150        ENDIF        ENDIF
151    
152        IF (selectMode .EQ. 0) THEN        IF (selectFindRoSurf.EQ.0 .AND. selectMode .GT. 0 ) THEN
153  C- Find Po_surf : Linear between 2 half levels :  C- Find Po_surf : Linear between 2 half levels :
154        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
155         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
# Line 168  c          IF (Hfld(i,j,bi,bj).GT.phiRef Line 171  c          IF (Hfld(i,j,bi,bj).GT.phiRef
171        ENDDO        ENDDO
172    
173  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
174  C-- endif selectMode=0  C-- endif selectFindRoSurf=0 & selectMode > 0
175        ENDIF        ENDIF
176    
177        IF (abs(selectMode) .EQ. 1) THEN        IF ( selectFindRoSurf.EQ.1 ) THEN
178  C-- define high resolution Pressure discretization:  C-- define high resolution Pressure discretization:
179    
180        recip_kappa = 1. _d 0 / atm_kappa        recip_kappa = 1. _d 0 / atm_kappa
# Line 187  C-- define high resolution Pressure disc Line 190  C-- define high resolution Pressure disc
190          dPiHvR(k) = PiHvR(k) - PiHvR(k+1)          dPiHvR(k) = PiHvR(k) - PiHvR(k+1)
191        ENDDO        ENDDO
192    
193    C--   to modify pressure when using non fully linear relation between Phi & p
194    C       (Integr_GeoPot=2 & Tracer Point at middle between 2 interfaces)
195          DO k=1,Nr
196             ratioRm(k) = 1. _d 0
197             ratioRp(k) = 1. _d 0
198             IF (k.GT.1 ) ratioRm(k) = 0.5 _d 0*drC(k)/(rF(k)-rC(k))
199             IF (k.LT.Nr) ratioRp(k) = 0.5 _d 0*drC(k+1)/(rC(k)-rF(k+1))
200          ENDDO
201    
202  #ifdef CHECK_ANALYLIC_THETA  #ifdef CHECK_ANALYLIC_THETA
203        _BEGIN_MASTER( mythid )        _BEGIN_MASTER( mythid )
204        DO j=1,61        DO j=1,61
# Line 202  C-- define high resolution Pressure disc Line 214  C-- define high resolution Pressure disc
214        STOP 'CHECK_ANALYLIC_THETA'        STOP 'CHECK_ANALYLIC_THETA'
215  #endif /* CHECK_ANALYLIC_THETA */  #endif /* CHECK_ANALYLIC_THETA */
216    
217  C-- endif |selectMode|=1  C-- endif selectFindRoSurf=1
218        ENDIF        ENDIF
219    
220        IF (selectMode .EQ. 1) THEN        IF ( selectFindRoSurf*selectMode .GT. 0) THEN
221  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
222  C- Find Po_surf such as g*Hfld = Phi[Po_surf,theta(yLat,p)]:  C- Find Po_surf such as g*Hfld = Phi[Po_surf,theta(yLat,p)]:
223    
224        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
225         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
226  C- start bi,bj loop:  C- start bi,bj loop:
227    
228          DO j=1,sNy          DO j=1,sNy
229           DO i=1,sNx           DO i=1,sNx
230            IF ( Hfld(i,j,bi,bj) .LE. 0. _d 0) THEN            IF ( Hfld(i,j,bi,bj) .LE. 0. _d 0) THEN
# Line 226  C- start bi,bj loop: Line 239  C- start bi,bj loop:
239  C-    compute  phi/g corresponding to next p_level:  C-    compute  phi/g corresponding to next p_level:
240               dzLoc = dPiHvR(k)*thetaHvR(k)*recip_gravity               dzLoc = dPiHvR(k)*thetaHvR(k)*recip_gravity
241               IF ( Hfld(i,j,bi,bj) .LE. zLoc+dzLoc ) THEN               IF ( Hfld(i,j,bi,bj) .LE. zLoc+dzLoc ) THEN
242  C-    compute the normalized surf. Pressure Po_surf  C-    compute the normalized surf. Pressure psNorm
243                 PiLoc = PiHvR(k)                 PiLoc = PiHvR(k)
244       &               - gravity*(Hfld(i,j,bi,bj)-zLoc)/thetaHvR(k)       &               - gravity*(Hfld(i,j,bi,bj)-zLoc)/thetaHvR(k)
245                 Po_surf = (PiLoc/atm_Cp)**recip_kappa                 psNorm = (PiLoc/atm_Cp)**recip_kappa
246  C- use linear interpolation:  C- use linear interpolation:
247  c              Po_surf = pLevHvR(k)  c              psNorm = pLevHvR(k)
248  c    &                 - dpHvR*(Hfld(i,j,bi,bj)-zLoc)/dzLoc  c    &                - dpHvR*(Hfld(i,j,bi,bj)-zLoc)/dzLoc
249                 zLoc = -1.                 zLoc = -1.
250               ELSE               ELSE
251                 zLoc = zLoc + dzLoc                 zLoc = zLoc + dzLoc
# Line 250  c    &                 - dpHvR*(Hfld(i,j Line 263  c    &                 - dpHvR*(Hfld(i,j
263               CALL PRINT_ERROR( msgBuf , myThid)               CALL PRINT_ERROR( msgBuf , myThid)
264               STOP 'ABNORMAL END: S/R INI_P_GROUND'                 STOP 'ABNORMAL END: S/R INI_P_GROUND'  
265             ELSE             ELSE
266               Pfld(i,j,bi,bj) = Po_surf*atm_Po               Pfld(i,j,bi,bj) = psNorm*atm_Po
267             ENDIF             ENDIF
268            ENDIF            ENDIF
269           ENDDO           ENDDO
270          ENDDO          ENDDO
271    
272            IF (selectMode.EQ.2 .AND. integr_GeoPot.NE.1) THEN
273    C---------
274    C     Modify pressure to account for non fully linear relation between
275    C      Phi & p (due to numerical truncation of the Finite Difference scheme)
276    C---------
277              DO j=1,sNy
278               DO i=1,sNx
279                 Po_surf = Pfld(i,j,bi,bj)
280                  IF ( Po_surf.LT.rC(1) .AND. Po_surf.GT.rC(Nr) ) THEN
281                    findPoSurf = .TRUE.
282                    DO k=1,Nr
283                      IF ( findPoSurf .AND. Po_surf.GE.rC(k) ) THEN
284                        Po_surf = rC(k) + (Po_surf-rC(k))/ratioRm(k)
285                        findPoSurf = .FALSE.
286                      ENDIF
287                      rMidKp1 = rF(k+1)
288                      IF (k.LT.Nr) rMidKp1 = (rC(k)+rC(k+1))*0.5 _d 0
289                      IF ( findPoSurf .AND. Po_surf.GE.rMidKp1 ) THEN
290                        Po_surf = rC(k) + (Po_surf-rC(k))/ratioRp(k)
291                        findPoSurf = .FALSE.
292                      ENDIF
293                    ENDDO
294                    IF ( findPoSurf )
295         &               STOP 'S/R INI_P_GROUND: Pb with selectMode=2'
296                  ENDIF
297                 Pfld(i,j,bi,bj) = Po_surf
298               ENDDO
299              ENDDO
300    C---------
301            ENDIF
302    
303  C- end bi,bj loop.  C- end bi,bj loop.
304         ENDDO         ENDDO
305        ENDDO        ENDDO
306    
307  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
308  C-- endif selectMode=1  C-- endif selectFindRoSurf*selectMode > 0
309        ENDIF        ENDIF
310    
311        IF (selectMode .EQ. -1) THEN        IF (selectMode .LT. 0) THEN
312  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
313  C  goal: evaluate phi0surf (used for computing geopotential'= Phi - PhiRef):  C--- Compute Hfld = Phi(Pfld)/g.
 C   phi0surf = g*Ztopo-PhiRef(Ro_surf) if no truncation was applied to Ro_surf;  
 C  but because of hFacMin, surf.ref.pressure (=Ro_surf) is truncated and  
 C   phi0surf = Phi(Theta-Analytic,P=Ro_surf) - PhiRef(P=Ro_surf)  
 C-----  
314    
315        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
316         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
317  C- start bi,bj loop:  C- start bi,bj loop:
318    
319  C--   Compute Hfld from g*Hfld = Phi[Po_surf,theta(yLat,p)]:  C--  Compute Hfld from g*Hfld = PhiRef(Po_surf)
         DO j=1,sNy  
          DO i=1,sNx  
           IF ( Pfld(i,j,bi,bj) .GE. Ro_SeaLevel) THEN  
            Hfld(i,j,bi,bj) = 0.  
           ELSE  
            Po_surf = Pfld(i,j,bi,bj)/atm_Po  
            kLev = 1 + INT( (pLevHvR(1)-Po_surf)/dpHvR )  
            yLatLoc  = yC(i,j,bi,bj)  
            CALL ANALYLIC_THETA( yLatLoc , pMidHvR,  
      &                       thetaHvR, kLev, mythid)  
 C-    compute height at level pLev(kLev) just below Pfld:  
            zLoc = 0.  
            DO k=1,kLev-1  
              dzLoc = dPiHvR(k)*thetaHvR(k)*recip_gravity  
              zLoc = zLoc + dzLoc  
            ENDDO  
            dzLoc = ( PiHvR(kLev)-atm_Cp*(Po_surf**atm_kappa) )  
      &         * thetaHvR(kLev)*recip_gravity  
            Hfld(i,j,bi,bj) = zLoc + dzLoc  
           ENDIF  
          ENDDO  
         ENDDO  
   
 C--  compute phi0surf (stored in Hfld):  
320          DO j=1,sNy          DO j=1,sNy
321           DO i=1,sNx           DO i=1,sNx
322  C-   compute phiLoc = PhiRef(Ro_surf):  C-   compute phiLoc = PhiRef(Ro_surf):
# Line 314  C-   compute phiLoc = PhiRef(Ro_surf): Line 331  C-   compute phiLoc = PhiRef(Ro_surf):
331       &       + (phiRef(2*ks+1)-phiRef(2*ks))       &       + (phiRef(2*ks+1)-phiRef(2*ks))
332       &        *(Pfld(i,j,bi,bj)-rC(ks))/(rHalf(2*ks+1)-rHalf(2*ks))       &        *(Pfld(i,j,bi,bj)-rC(ks))/(rHalf(2*ks+1)-rHalf(2*ks))
333             ENDIF             ENDIF
334             Hfld(i,j,bi,bj) = gravity*(Hfld(i,j,bi,bj) - phiLoc)             Hfld(i,j,bi,bj) = phiLoc
335            ELSE            ELSE
336             Hfld(i,j,bi,bj) = 0.             Hfld(i,j,bi,bj) = 0.
337            ENDIF            ENDIF
338           ENDDO           ENDDO
339          ENDDO          ENDDO
340    
341            IF (selectFindRoSurf.EQ.1) THEN
342    C-----
343    C  goal: evaluate phi0surf (used for computing geopotential'= Phi - PhiRef):
344    C   phi0surf = g*Ztopo-PhiRef(Ro_surf) if no truncation was applied to Ro_surf;
345    C  but because of hFacMin, surf.ref.pressure (=Ro_surf) is truncated and
346    C   phi0surf = Phi(Theta-Analytic,P=Ro_surf) - PhiRef(P=Ro_surf)
347    C-----
348    C--   Compute Hfld from g*Hfld = Phi[Po_surf,theta(yLat,p)]:
349             DO j=1,sNy
350              DO i=1,sNx
351               zLoc = 0.
352               IF ( Pfld(i,j,bi,bj) .LT. Ro_SeaLevel) THEN
353                Po_surf = Pfld(i,j,bi,bj)
354    C---------
355    C     Modify pressure to account for non fully linear relation between
356    C      Phi & p (due to numerical truncation of the Finite Difference scheme)
357                 IF (selectMode.EQ.-2 .AND. integr_GeoPot.NE.1) THEN
358                  IF ( Po_surf.LT.rC(1) .AND. Po_surf.GT.rC(Nr) ) THEN
359                    findPoSurf = .TRUE.
360                    DO k=1,Nr
361                      IF ( findPoSurf .AND. Po_surf.GE.rC(k) ) THEN
362                        Po_surf = rC(k) + (Po_surf-rC(k))*ratioRm(k)
363                        findPoSurf = .FALSE.
364                      ENDIF
365                      IF ( findPoSurf .AND. Po_surf.GE.rF(k+1) ) THEN
366                        Po_surf = rC(k) + (Po_surf-rC(k))*ratioRp(k)
367                        findPoSurf = .FALSE.
368                      ENDIF
369                    ENDDO
370                  ENDIF
371                 ENDIF
372    C---------
373                psNorm = Po_surf/atm_Po
374                kLev = 1 + INT( (pLevHvR(1)-psNorm)/dpHvR )
375                yLatLoc  = yC(i,j,bi,bj)
376                CALL ANALYLIC_THETA( yLatLoc , pMidHvR,
377         &                        thetaHvR, kLev, mythid)
378    C-    compute height at level pLev(kLev) just below Pfld:
379                DO k=1,kLev-1
380                  dzLoc = dPiHvR(k)*thetaHvR(k)*recip_gravity
381                  zLoc = zLoc + dzLoc
382                ENDDO
383                dzLoc = ( PiHvR(kLev)-atm_Cp*(psNorm**atm_kappa) )
384         &            * thetaHvR(kLev)*recip_gravity
385                zLoc = zLoc + dzLoc
386               ENDIF
387    C-    compute phi0surf = Phi[Po_surf,theta(yLat,p)] - PhiRef(Po_surf)
388               phi0surf(i,j,bi,bj) = gravity*(zLoc - Hfld(i,j,bi,bj))
389    C-    save Phi[Po_surf,theta(yLat,p)] in Hfld (replacing PhiRef(Po_surf)):
390               Hfld(i,j,bi,bj) = zLoc
391              ENDDO
392             ENDDO
393    C- endif selectFindRoSurf=1
394            ENDIF
395    
396  C- end bi,bj loop.  C- end bi,bj loop.
397         ENDDO         ENDDO
398        ENDDO        ENDDO
399    
400  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
401  C-- endif selectMode=-1  C-- endif selectMode < 0
402        ENDIF        ENDIF
403    
404        RETURN        RETURN
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.5
  ViewVC Help
Powered by ViewVC 1.1.22