/[MITgcm]/MITgcm_contrib/dgoldberg/streamice/streamice_advect_thickness.F
ViewVC logotype

Diff of /MITgcm_contrib/dgoldberg/streamice/streamice_advect_thickness.F

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

revision 1.1 by heimbach, Thu Mar 29 15:59:21 2012 UTC revision 1.3 by dgoldberg, Thu Jul 26 16:13:18 2012 UTC
# Line 24  C     === Global variables === Line 24  C     === Global variables ===
24  #include "PARAMS.h"  #include "PARAMS.h"
25  #include "STREAMICE.h"  #include "STREAMICE.h"
26  #include "STREAMICE_ADV.h"  #include "STREAMICE_ADV.h"
27    #ifdef ALLOW_AUTODIFF_TAMC
28    # include "tamc.h"
29    #endif
30    
31        INTEGER myThid        INTEGER myThid
32        _RL time_step        _RL time_step
# Line 34  C     === Global variables === Line 37  C     === Global variables ===
37        _RL thick_bd        _RL thick_bd
38        _RL SLOPE_LIMITER        _RL SLOPE_LIMITER
39        _RL sec_per_year, time_step_loc        _RL sec_per_year, time_step_loc
40          CHARACTER*(MAX_LEN_MBUF) msgBuf
41        external SLOPE_LIMITER        external SLOPE_LIMITER
42    
43        sec_per_year = 365.*86400.        sec_per_year = 365.*86400.
44    
45        time_step_loc = time_step / sec_per_year        time_step_loc = time_step / sec_per_year
46    
47    #ifdef ALLOW_AUTODIFF_TAMC
48    CADJ STORE streamice_hmask  = comlev1, key=ikey_dynamics
49    #endif
50    
51        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
52         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
53          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
# Line 62  C     === Global variables === Line 70  C     === Global variables ===
70         ENDDO         ENDDO
71        ENDDO        ENDDO
72    
73  !       PRINT *, "H in last row ", H_streamice(81,20,1,1)  
74    
75    #ifdef ALLOW_AUTODIFF_TAMC
76    CADJ STORE h_after_uflux_si = comlev1, key=ikey_dynamics
77    CADJ STORE streamice_hmask  = comlev1, key=ikey_dynamics
78    #endif
79    
80        CALL STREAMICE_ADVECT_THICKNESS_X ( myThid,        CALL STREAMICE_ADVECT_THICKNESS_X ( myThid,
81       O   hflux_x_SI,       O   hflux_x_SI,
82       O   h_after_uflux_SI,       O   h_after_uflux_SI,
83       I   time_step_loc )       I   time_step_loc )
84    
85  !       PRINT *, "H in last row ", H_streamice(81,20,1,1)  
86    
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)
# Line 82  C     === Global variables === Line 95  C     === Global variables ===
95         ENDDO         ENDDO
96        ENDDO        ENDDO
97    
98    #ifdef ALLOW_AUTODIFF_TAMC
99    CADJ STORE h_after_vflux_si = comlev1, key=ikey_dynamics
100    CADJ STORE streamice_hmask  = comlev1, key=ikey_dynamics
101    #endif
102    
103        CALL STREAMICE_ADVECT_THICKNESS_Y ( myThid,        CALL STREAMICE_ADVECT_THICKNESS_Y ( myThid,
104       O   hflux_y_SI,       O   hflux_y_SI,
105       O   h_after_vflux_SI,       O   h_after_vflux_SI,
106       I   time_step_loc )       I   time_step_loc )
107    
108  !       PRINT *, "H in last row ", H_streamice(81,20,1,1)  
109    
110        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
111         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
# Line 102  C     === Global variables === Line 120  C     === Global variables ===
120         ENDDO         ENDDO
121        ENDDO        ENDDO
122    
123  !       PRINT *, "H in last row ", H_streamice(81,20,1,1)  
124    
125        CALL STREAMICE_ADV_FRONT ( myThid, time_step_loc )        CALL STREAMICE_ADV_FRONT ( myThid, time_step_loc )
126    
127  !       PRINT *, "H in last row ", H_streamice(81,20,1,1)  
128    
129        _EXCH_XY_RL( H_streamice, myThid )        _EXCH_XY_RL( H_streamice, myThid )
130        _EXCH_XY_RL( area_shelf_streamice, myThid )        _EXCH_XY_RL( area_shelf_streamice, myThid )
131        _EXCH_XY_RL( STREAMICE_hmask, myThid )        _EXCH_XY_RL( STREAMICE_hmask, myThid )
132                
133        PRINT *, "END STREAMICE_ADVECT_THICKNESS"        WRITE(msgBuf,'(A)') 'END STREAMICE_ADVECT_THICKNESS'
134           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
135  #endif       &                     SQUEEZE_RIGHT , 1)
       RETURN  
       END SUBROUTINE STREAMICE_ADVECT_THICKNESS  
   
 !     NEED TO ADD SOME SORT OF CHECK THAT THE HALOS ARE LARGE ENOUGH  
   
       SUBROUTINE STREAMICE_ADVECT_THICKNESS_X ( myThid ,      
      O   hflux_x ,  
      O   h ,  
      I   time_step )  
   
       IMPLICIT NONE  
   
 C     O   hflux_x ! flux per unit width across face  
 C     O   h  
 C     I   time_step  
   
 C     === Global variables ===  
 #include "SIZE.h"  
 #include "GRID.h"  
 #include "EEPARAMS.h"  
 #include "PARAMS.h"  
 #include "STREAMICE.h"  
   
       INTEGER myThid  
       _RL hflux_x          (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)  
       _RL h                (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)    
       _RL time_step  
136                
 #ifdef ALLOW_STREAMICE  
   
 C     LOCAL VARIABLES  
   
       INTEGER i, j, bi, bj, Gi, Gj, k  
       _RL uface, phi  
       _RL stencil (-1:1)  
       LOGICAL H0_valid  ! there are valid cells to calculate a  
                         ! slope-limited 2nd order flux  
       _RL SLOPE_LIMITER  
       _RL total_vol_out  
       external SLOPE_LIMITER  
   
        total_vol_out = 0.0  
   
       DO bj=myByLo(myThid),myByHi(myThid)  
        DO bi=myBxLo(myThid),myBxHi(myThid)  
         DO j=1-3,sNy+3  
          Gj = (myYGlobalLo-1)+(bj-1)*sNy+j  
          IF ((Gj .ge. 1) .and. (Gj .le. Ny)) THEN  
           DO i=1-2,sNx+3  
 C        THESE ARRAY BOUNDS INSURE THAT AFTER THIS STEP,  
 C        VALUES WILL BE RELIABLE 2 GRID CELLS OUT IN THE  
 C        X DIRECTION AND 3 CELLS OUT IN THE Y DIR  
            IF ((STREAMICE_hmask(i,j,bi,bj).eq.1.0) .or.  
      &         ((STREAMICE_hmask(i-1,j,bi,bj).eq.1.0) .and.  
      &          (STREAMICE_hmask(i,j,bi,bj).ne.1.0))) THEN  
   
             Gi = (myXGlobalLo-1)+(bi-1)*sNx+i  
               
             IF (STREAMICE_ufacemask(i,j,bi,bj).eq.4.0) THEN  
              hflux_x (i,j,bi,bj) = u_flux_bdry_SI (i,j,bi,bj)  
             ELSE  
             
              uface = .5 *  
      &        (U_streamice(i,j,bi,bj)+U_streamice(i,j+1,bi,bj))  
   
   
              IF (uface .gt. 0. _d 0) THEN  
               DO k=-1,1  
                stencil (k) = h(i+k-1,j,bi,bj)  
               ENDDO  
               IF ((STREAMICE_hmask(i,j,bi,bj).eq.1.0) .and.  
      &            (STREAMICE_hmask(i-2,j,bi,bj).eq.1.0)) H0_valid=.true.  
               
               IF ((Gi.eq.1).and.(STREAMICE_hmask(i-1,j,bi,bj).eq.3.0))  
      &         THEN  ! we are at western bdry and there is a thick. bdry cond  
                hflux_x (i,j,bi,bj) = h(i-1,j,bi,bj) * uface  
               ELSEIF (H0_valid) THEN  
                phi = SLOPE_LIMITER (  
      &          stencil(0)-stencil(-1),  
      &          stencil(1)-stencil(0))  
                hflux_x (i,j,bi,bj) = uface *  
      &          (stencil(0) - phi * .5 * (stencil(0)-stencil(1)))  
               ELSE ! one of the two cells needed for a HO scheme is missing, use FO scheme  
                hflux_x (i,j,bi,bj) = uface * stencil(0)  
               ENDIF  
               
              ELSE ! uface <= 0  
   
               DO k=-1,1  
                stencil (k) = h(i-k,j,bi,bj)  
               ENDDO  
               IF ((STREAMICE_hmask(i-1,j,bi,bj).eq.1.0) .and.  
      &            (STREAMICE_hmask(i+1,j,bi,bj).eq.1.0)) H0_valid=.true.  
   
               IF ((Gi.eq.Nx).and.(STREAMICE_hmask(i+1,j,bi,bj).eq.3.0))  
      &         THEN  ! we are at western bdry and there is a thick. bdry cond  
                hflux_x (i,j,bi,bj) = h(i+1,j,bi,bj) * uface  
               ELSEIF (H0_valid) THEN  
                phi = SLOPE_LIMITER (  
      &          stencil(0)-stencil(-1),  
      &          stencil(1)-stencil(0))  
                hflux_x (i,j,bi,bj) = uface *  
      &          (stencil(0) - phi * .5 * (stencil(0)-stencil(1)))  
               ELSE ! one of the two cells needed for a HO scheme is missing, use FO scheme  
                hflux_x (i,j,bi,bj) = uface * stencil(0)  
               ENDIF  
   
              ENDIF  
   
             ENDIF  
   
             if (streamice_ufacemask(i,j,bi,bj).eq.2.0) THEN  
              total_vol_out = total_vol_out +  
      &         hflux_x (i,j,bi,bj) * time_step  
             ENDIF  
   
            ENDIF  
           ENDDO  
          ENDIF  
         ENDDO  
        ENDDO  
       ENDDO  
   
 C     X-FLUXES AT CELL BOUNDARIES CALCULATED; NOW TAKE FLUX DIVERGENCE TO INCREMENT THICKNESS  
         
       DO bj=myByLo(myThid),myByHi(myThid)  
        DO bi=myBxLo(myThid),myBxHi(myThid)  
         DO j=1-3,sNy+3  
          Gj = (myYGlobalLo-1)+(bj-1)*sNy+j  
          IF ((Gj .ge. 1) .and. (Gj .le. Ny)) THEN  
           DO i=1-2,sNx+2  
            IF (STREAMICE_hmask(i,j,bi,bj).eq.1.0) THEN  
             h(i,j,bi,bj) = h(i,j,bi,bj) - time_step *  
      &       (hflux_x(i+1,j,bi,bj)*dyG(i+1,j,bi,bj) -  
      &        hflux_x(i,j,bi,bj)*dyG(i,j,bi,bj)) *  
      &       recip_rA (i,j,bi,bj)  
            ENDIF  
           ENDDO  
          ENDIF  
         ENDDO  
        ENDDO  
       ENDDO  
   
 !       PRINT *, "TOTAL VOLUME OUT: ", total_vol_out  
         
 #endif  
       RETURN  
       END SUBROUTINE STREAMICE_ADVECT_THICKNESS_X  
   
       SUBROUTINE STREAMICE_ADVECT_THICKNESS_Y ( myThid ,      
      O   hflux_y ,  
      O   h ,  
      I   time_step )  
   
       IMPLICIT NONE  
   
 C     O   hflux_y ! flux per unit width across face  
 C     O   h  
 C     I   time_step  
   
 C     === Global variables ===  
 #include "SIZE.h"  
 #include "GRID.h"  
 #include "EEPARAMS.h"  
 #include "PARAMS.h"  
 #include "STREAMICE.h"  
   
       INTEGER myThid  
       _RL hflux_y          (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)  
       _RL h                (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)    
       _RL time_step  
         
 #ifdef ALLOW_STREAMICE  
   
 C     LOCAL VARIABLES  
   
       INTEGER i, j, bi, bj, Gi, Gj, k  
       _RL vface, phi  
       _RL stencil (-1:1)  
       LOGICAL H0_valid  ! there are valid cells to calculate a  
                         ! slope-limited 2nd order flux  
       _RL SLOPE_LIMITER  
       external SLOPE_LIMITER  
   
       DO bj=myByLo(myThid),myByHi(myThid)  
        DO bi=myBxLo(myThid),myBxHi(myThid)  
         DO j=1-1,sNy+2  
          Gj = (myYGlobalLo-1)+(bj-1)*sNy+j  
          DO i=1-2,sNx+2  
           Gi = (myXGlobalLo-1)+(bi-1)*sNx+i  
           IF ((Gi.ge.1) .and. (Gi.le.Nx)) THEN  
 C         THESE ARRAY BOUNDS INSURE THAT AFTER THIS STEP,  
 C         VALUES WILL BE RELIABLE 1 GRID CELLS OUT IN THE  
 C         Y DIRECTION  
            IF ((STREAMICE_hmask(i,j,bi,bj).eq.1.0) .or.  
      &         ((STREAMICE_hmask(i,j-1,bi,bj).eq.1.0) .and.  
      &          (STREAMICE_hmask(i,j,bi,bj).ne.1.0))) THEN  
   
             IF (STREAMICE_vfacemask(i,j,bi,bj).eq.4.0) THEN  
              hflux_y (i,j,bi,bj) = v_flux_bdry_SI (i,j,bi,bj)  
             ELSE  
             
              vface = .5 *  
      &        (V_streamice(i,j,bi,bj)+V_streamice(i+1,j,bi,bj))  
   
              IF (vface .gt. 0. _d 0) THEN  
               DO k=-1,1  
                stencil (k) = h(i,j+k-1,bi,bj)  
               ENDDO  
               IF ((STREAMICE_hmask(i,j,bi,bj).eq.1.0) .and.  
      &            (STREAMICE_hmask(i,j-2,bi,bj).eq.1.0)) H0_valid=.true.  
               
               IF ((Gj.eq.1).and.(STREAMICE_hmask(i,j-1,bi,bj).eq.3.0))  
      &         THEN  ! we are at western bdry and there is a thick. bdry cond  
                hflux_y (i,j,bi,bj) = h(i,j-1,bi,bj) * vface  
               ELSEIF (H0_valid) THEN  
                phi = SLOPE_LIMITER (  
      &          stencil(0)-stencil(-1),  
      &          stencil(1)-stencil(0))  
                hflux_y (i,j,bi,bj) = vface *  
      &          (stencil(0) - phi * .5 * (stencil(0)-stencil(1)))  
               ELSE ! one of the two cells needed for a HO scheme is missing, use FO scheme  
                hflux_y (i,j,bi,bj) = vface * stencil(0)  
               ENDIF  
              ELSE ! uface <= 0  
               DO k=-1,1  
                stencil (k) = h(i,j-k,bi,bj)  
               ENDDO  
               IF ((STREAMICE_hmask(i,j-1,bi,bj).eq.1.0) .and.  
      &            (STREAMICE_hmask(i,j+1,bi,bj).eq.1.0)) H0_valid=.true.  
   
               IF ((Gj.eq.Ny).and.(STREAMICE_hmask(i,j+1,bi,bj).eq.3.0))  
      &         THEN  ! we are at western bdry and there is a thick. bdry cond  
                hflux_y (i,j,bi,bj) = h(i,j+1,bi,bj) * vface  
               ELSEIF (H0_valid) THEN  
                phi = SLOPE_LIMITER (  
      &          stencil(0)-stencil(-1),  
      &          stencil(1)-stencil(0))  
                hflux_y (i,j,bi,bj) = vface *  
      &          (stencil(0) - phi * .5 * (stencil(0)-stencil(1)))  
               ELSE ! one of the two cells needed for a HO scheme is missing, use FO scheme  
                hflux_y (i,j,bi,bj) = vface * stencil(0)  
               ENDIF  
   
              ENDIF ! uface  0  
   
             ENDIF  
            ENDIF  
           ENDIF  
          ENDDO  
         ENDDO  
        ENDDO  
       ENDDO  
   
 C     X-FLUXES AT CELL BOUNDARIES CALCULATED; NOW TAKE FLUX DIVERGENCE TO INCREMENT THICKNESS  
   
         
         
       DO bj=myByLo(myThid),myByHi(myThid)  
        DO bi=myBxLo(myThid),myBxHi(myThid)  
         DO j=1-1,sNy+1  
          DO i=1-2,sNx+2  
           Gi = (myXGlobalLo-1)+(bi-1)*sNx+i  
           IF ((Gi .ge. 1) .and. (Gi .le. Nx)) THEN  
            IF (STREAMICE_hmask(i,j,bi,bj).eq.1.0) THEN  
             h(i,j,bi,bj) = h(i,j,bi,bj) - time_step *  
      &       (hflux_y(i,j+1,bi,bj)*dxG(i,j+1,bi,bj) -  
      &        hflux_y(i,j,bi,bj)*dxG(i,j,bi,bj)) *  
      &       recip_rA (i,j,bi,bj)  
            ENDIF  
           ENDIF  
          ENDDO  
         ENDDO  
        ENDDO  
       ENDDO  
   
 !       CALL WRITE_FLD_XY_RL ("h_after_yflux","",  
 !      & h, 0, myThid)  
   
137  #endif  #endif
138        RETURN        END
       END SUBROUTINE STREAMICE_ADVECT_THICKNESS_Y  
139    
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.3
  ViewVC Help
Powered by ViewVC 1.1.22