Back to home page

MITgcm
 
 

    

File indexing completed on 2018-03-02 18:44:16 UTC

view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
5ca83cd8f7 Dani*0001 #include "STREAMICE_OPTIONS.h"
                0002 
                0003 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0004 
                0005 CBOP
                0006       SUBROUTINE STREAMICE_DRIVING_STRESS_PPM( myThid ) 
                0007 !      O taudx, 
                0008 !      O taudy )
                0009 
                0010 C     /============================================================\
                0011 C     | SUBROUTINE                                                 |   
                0012 C     | o                                                          |
                0013 C     |============================================================|
                0014 C     |                                                            |
                0015 C     \============================================================/
                0016       IMPLICIT NONE
                0017 
                0018 C     === Global variables ===
                0019 #include "SIZE.h"
                0020 #include "EEPARAMS.h"
                0021 #include "PARAMS.h"
                0022 #include "GRID.h"
                0023 #include "STREAMICE.h"
                0024 #include "STREAMICE_CG.h"
                0025 
                0026 C     !INPUT/OUTPUT ARGUMENTS
                0027       INTEGER myThid
                0028 !       _RL taudx (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0029 !       _RL taudx (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0030 
                0031 #ifdef ALLOW_STREAMICE
                0032 
                0033 
                0034 C     LOCAL VARIABLES
                0035       INTEGER i, j, bi, bj, k, l,
                0036      &        Gi, Gj
                0037       LOGICAL at_west_bdry, at_east_bdry, 
                0038      &        at_north_bdry, at_south_bdry
                0039       _RL sx, sy, diffx, diffy, neu_val
                0040       
                0041       IF (myXGlobalLo.eq.1) at_west_bdry = .true.
                0042       IF (myYGlobalLo.eq.1) at_south_bdry = .true.
                0043       IF (myXGlobalLo-1+sNx*nSx.eq.Nx) 
                0044      & at_east_bdry = .false.
                0045       IF (myYGlobalLo-1+sNy*nSy.eq.Ny) 
                0046      & at_north_bdry = .false.
                0047 
                0048       DO bj = myByLo(myThid), myByHi(myThid)
                0049        DO bi = myBxLo(myThid), myBxHi(myThid)
                0050         DO j=1-OLy,sNy+OLy
                0051          DO i=1-OLx,sNx+OLx
                0052           taudx_SI(i,j,bi,bj) = 0. _d 0
                0053           taudy_SI(i,j,bi,bj) = 0. _d 0
                0054          ENDDO
                0055         ENDDO
                0056        ENDDO
                0057       ENDDO
                0058 
                0059       DO bj = myByLo(myThid), myByHi(myThid)
                0060        DO bi = myBxLo(myThid), myBxHi(myThid)
                0061         
                0062         DO i=0,sNx+1
                0063          DO j=0,sNy+1
                0064 
                0065           diffx = 0. _d 0
                0066           diffy = 0. _d 0
                0067           sx = 0. _d 0
                0068           sy = 0. _d 0
                0069 
                0070           Gi = (myXGlobalLo-1)+(bi-1)*sNx+i
                0071           Gj = (myYGlobalLo-1)+(bj-1)*sNy+j
                0072 
                0073           IF (STREAMICE_hmask(i,j,bi,bj).eq.1.0) THEN
                0074 
                0075            ! we are in an "active" cell
                0076 
                0077            IF (Gi.eq.1.AND..NOT.STREAMICE_EW_periodic) THEN
                0078 
                0079             ! western boundary - only one sided possible
                0080 
                0081             IF (STREAMICE_hmask(i+1,j,bi,bj).eq.1.0) THEN
                0082 
                0083              ! cell to east is active
                0084 
                0085              sx = (surf_el_streamice(i+1,j,bi,bj)-
                0086      &             surf_el_streamice(i,j,bi,bj))/dxC(i+1,j,bi,bj)
                0087             ELSE
                0088 
                0089              ! cell to east is empty
                0090 
                0091              sx = 0. _d 0
                0092             ENDIF
                0093 
                0094             DO k=0,1
                0095              DO l=0,1
                0096               IF (STREAMICE_umask(i+k,j+l,bi,bj).eq.1.0) THEN
                0097                taudx_SI(i+k,j+l,bi,bj) = taudx_SI(i+k,j+l,bi,bj) -
                0098      &          0.25 * streamice_density * gravity * 
                0099      &          (streamice_bg_surf_slope_x+sx) * 
                0100      &          H_streamice(i,j,bi,bj) * rA(i,j,bi,bj)
                0101               ENDIF
                0102              ENDDO
                0103             ENDDO
                0104 
                0105            ELSEIF (Gi.eq.Nx.AND..NOT.STREAMICE_EW_periodic) THEN
                0106 
                0107             ! eastern boundary - only one sided possible
                0108 
                0109             IF (STREAMICE_hmask(i-1,j,bi,bj).eq.1.0) THEN
                0110 
                0111              ! cell to west is active
                0112 
                0113              sx = (surf_el_streamice(i,j,bi,bj)-
                0114      &             surf_el_streamice(i-1,j,bi,bj))/dxC(i,j,bi,bj)
                0115             ELSE
                0116 
                0117              ! cell to west is inactive
                0118 
                0119              sx = 0. _d 0
                0120             ENDIF
                0121 
                0122             DO k=0,1
                0123              DO l=0,1
                0124               IF (STREAMICE_umask(i+k,j+l,bi,bj).eq.1.0) THEN
                0125                taudx_SI(i+k,j+l,bi,bj) = taudx_SI(i+k,j+l,bi,bj) -
                0126      &          0.25 * streamice_density * gravity * 
                0127      &          (streamice_bg_surf_slope_x+sx) * 
                0128      &          H_streamice(i,j,bi,bj) * rA(i,j,bi,bj)
                0129               ENDIF
                0130              ENDDO
                0131             ENDDO
                0132 
                0133            ELSE
                0134 
                0135             ! interior (west-east) cell
                0136 
                0137             IF (STREAMICE_hmask(i+1,j,bi,bj).eq.1.0 .and.
                0138      &          STREAMICE_hmask(i-1,j,bi,bj).eq.1.0) THEN
                0139 
                0140              k = 0
                0141              DO l=0,1
                0142               IF (STREAMICE_umask(i+k,j+l,bi,bj).eq.1.0) THEN
                0143                taudx_SI(i+k,j+l,bi,bj) = taudx_SI(i+k,j+l,bi,bj) -
                0144      &          streamice_density * gravity * (1./6.) *
                0145      &          (-2.*surf_el_streamice(i-1,j,bi,bj) + 
                0146      &           surf_el_streamice(i,j,bi,bj) + 
                0147      &           surf_el_streamice(i+1,j,bi,bj) +
                0148      &           3.*streamice_bg_surf_slope_x * dxF(i,j,bi,bj)) *
                0149      &          H_streamice(i,j,bi,bj) * .5 * dyF(i,j,bi,bj)
                0150               ENDIF
                0151              ENDDO
                0152 
                0153              k = 1
                0154              DO l=0,1
                0155               IF (STREAMICE_umask(i+k,j+l,bi,bj).eq.1.0) THEN
                0156                taudx_SI(i+k,j+l,bi,bj) = taudx_SI(i+k,j+l,bi,bj) -
                0157      &          streamice_density * gravity * (1./6.) *
                0158      &          (-surf_el_streamice(i-1,j,bi,bj) - 
                0159      &           surf_el_streamice(i,j,bi,bj) + 
                0160      &           2*surf_el_streamice(i+1,j,bi,bj) +
                0161      &           3.*streamice_bg_surf_slope_x * dxF(i,j,bi,bj)) *
                0162      &          H_streamice(i,j,bi,bj) * .5 * dyF(i,j,bi,bj)
                0163               ENDIF
                0164              ENDDO
                0165 
                0166 
                0167             ELSE 
                0168 
                0169              IF (STREAMICE_hmask(i+1,j,bi,bj).eq.1.0) THEN
                0170 
                0171               sx = (surf_el_streamice(i+1,j,bi,bj)-
                0172      &             surf_el_streamice(i,j,bi,bj))/dxC(i+1,j,bi,bj)
                0173 
                0174              ELSEIF (STREAMICE_hmask(i-1,j,bi,bj).eq.1.0) THEN
                0175 
                0176               sx = (surf_el_streamice(i,j,bi,bj)-
                0177      &             surf_el_streamice(i-1,j,bi,bj))/dxC(i,j,bi,bj)
                0178 
                0179              ELSE
                0180 
                0181               sx = 0. _d 0
                0182 
                0183              ENDIF
                0184 
                0185              DO k=0,1
                0186               DO l=0,1
                0187                IF (STREAMICE_umask(i+k,j+l,bi,bj).eq.1.0) THEN
                0188                 taudx_SI(i+k,j+l,bi,bj) = taudx_SI(i+k,j+l,bi,bj) -
                0189      &           0.25 * streamice_density * gravity * 
                0190      &           (streamice_bg_surf_slope_x+sx) * 
                0191      &           H_streamice(i,j,bi,bj) * rA(i,j,bi,bj)
                0192                ENDIF
                0193               ENDDO
                0194              ENDDO
                0195 
                0196             ENDIF
                0197 
                0198            ENDIF
                0199 
                0200 !!!!!!!! DONE WITH X-GRADIENT
                0201 
                0202            IF (Gj.eq.1.AND..NOT.STREAMICE_NS_periodic) THEN
                0203 
                0204             ! western boundary - only one sided possible
                0205 
                0206             IF (STREAMICE_hmask(i,j+1,bi,bj).eq.1.0) THEN
                0207 
                0208              ! cell to east is active
                0209 
                0210              sy = (surf_el_streamice(i,j+1,bi,bj)-
                0211      &             surf_el_streamice(i,j,bi,bj))/dyC(i,j+1,bi,bj)
                0212             ELSE
                0213 
                0214              ! cell to east is empty
                0215 
                0216              sy = 0. _d 0
                0217             ENDIF
                0218 
                0219             DO k=0,1
                0220              DO l=0,1
                0221               IF (STREAMICE_umask(i+k,j+l,bi,bj).eq.1.0) THEN
                0222                taudy_SI(i+k,j+l,bi,bj) = taudy_SI(i+k,j+l,bi,bj) -
                0223      &          0.25 * streamice_density * gravity * 
                0224      &          (streamice_bg_surf_slope_y+sy) * 
                0225      &          H_streamice(i,j,bi,bj) * rA(i,j,bi,bj)
                0226               ENDIF
                0227              ENDDO
                0228             ENDDO
                0229 
                0230            ELSEIF (Gj.eq.Ny.AND..NOT.STREAMICE_NS_periodic) THEN
                0231 
                0232             ! eastern boundary - only one sided possible
                0233 
                0234             IF (STREAMICE_hmask(i,j-1,bi,bj).eq.1.0) THEN
                0235 
                0236              ! cell to west is active
                0237 
                0238              sy = (surf_el_streamice(i,j,bi,bj)-
                0239      &             surf_el_streamice(i,j-1,bi,bj))/dyC(i,j,bi,bj)
                0240 
                0241             ELSE
                0242 
                0243              ! cell to west is inactive
                0244 
                0245              sy = 0. _d 0
                0246             ENDIF
                0247 
                0248             DO k=0,1
                0249              DO l=0,1
                0250               IF (STREAMICE_umask(i+k,j+l,bi,bj).eq.1.0) THEN
                0251                taudy_SI(i+k,j+l,bi,bj) = taudy_SI(i+k,j+l,bi,bj) -
                0252      &          0.25 * streamice_density * gravity * 
                0253      &          (streamice_bg_surf_slope_y+sy) * 
                0254      &          H_streamice(i,j,bi,bj) * rA(i,j,bi,bj)
                0255               ENDIF
                0256              ENDDO
                0257             ENDDO
                0258 
                0259            ELSE
                0260 
                0261             ! interior (west-east) cell
                0262 
                0263             IF (STREAMICE_hmask(i,j+1,bi,bj).eq.1.0 .and.
                0264      &          STREAMICE_hmask(i,j-1,bi,bj).eq.1.0) THEN
                0265 
                0266              l = 0
                0267              DO k=0,1
                0268               IF (STREAMICE_umask(i+k,j+l,bi,bj).eq.1.0) THEN
                0269                taudy_SI(i+k,j+l,bi,bj) = taudy_SI(i+k,j+l,bi,bj) -
                0270      &          streamice_density * gravity * (1./6.) *
                0271      &          (-2.*surf_el_streamice(i,j-1,bi,bj) + 
                0272      &           surf_el_streamice(i,j,bi,bj) + 
                0273      &           surf_el_streamice(i,j+1,bi,bj) +
                0274      &           3.*streamice_bg_surf_slope_y * dyF(i,j,bi,bj)) *
                0275      &          H_streamice(i,j,bi,bj) * .5 * dxF(i,j,bi,bj)
                0276               ENDIF
                0277              ENDDO
                0278 
                0279              l = 1
                0280              DO k=0,1
                0281               IF (STREAMICE_umask(i+k,j+l,bi,bj).eq.1.0) THEN
                0282                taudy_SI(i+k,j+l,bi,bj) = taudy_SI(i+k,j+l,bi,bj) -
                0283      &          streamice_density * gravity * (1./6.) *
                0284      &          (-surf_el_streamice(i,j-1,bi,bj) - 
                0285      &           surf_el_streamice(i,j,bi,bj) + 
                0286      &           2*surf_el_streamice(i,j+1,bi,bj) +
                0287      &           3.*streamice_bg_surf_slope_y * dyF(i,j,bi,bj)) *
                0288      &          H_streamice(i,j,bi,bj) * .5 * dxF(i,j,bi,bj)
                0289               ENDIF
                0290              ENDDO
                0291 
                0292 
                0293             ELSE 
                0294 
                0295              IF (STREAMICE_hmask(i,j+1,bi,bj).eq.1.0) THEN
                0296 
                0297               sy = (surf_el_streamice(i,j+1,bi,bj)-
                0298      &             surf_el_streamice(i,j,bi,bj))/dxC(i,j+1,bi,bj)
                0299 
                0300              ELSEIF (STREAMICE_hmask(i,j-1,bi,bj).eq.1.0) THEN
                0301 
                0302               sy = (surf_el_streamice(i,j,bi,bj)-
                0303      &             surf_el_streamice(i,j-1,bi,bj))/dxC(i,j,bi,bj)
                0304 
                0305              ELSE
                0306 
                0307               sy = 0. _d 0
                0308 
                0309              ENDIF
                0310 
                0311              DO k=0,1
                0312               DO l=0,1
                0313                IF (STREAMICE_umask(i+k,j+l,bi,bj).eq.1.0) THEN
                0314                 taudy_SI(i+k,j+l,bi,bj) = taudy_SI(i+k,j+l,bi,bj) -
                0315      &           0.25 * streamice_density * gravity * 
                0316      &           (streamice_bg_surf_slope_y+sy) * 
                0317      &           H_streamice(i,j,bi,bj) * rA(i,j,bi,bj)
                0318                ENDIF
                0319               ENDDO
                0320              ENDDO
                0321 
                0322             ENDIF
                0323 
                0324            ENDIF
                0325 
                0326 !            DO k=0,1
                0327 !             DO l=0,1
                0328 !              IF (STREAMICE_umask(i+k,j+l,bi,bj).eq.1.0) THEN
                0329 !               taudx_SI(i+k,j+l,bi,bj) = taudx_SI(i+k,j+l,bi,bj) -
                0330 !      &         0.25 * streamice_density * gravity * 
                0331 !      &         (streamice_bg_surf_slope_x+sx) * 
                0332 !      &         H_streamice(i,j,bi,bj) * rA(i,j,bi,bj)
                0333 !               taudy_SI(i+k,j+l,bi,bj) = taudy_SI(i+k,j+l,bi,bj) -
                0334 !      &         0.25 * streamice_density * gravity * 
                0335 !      &         (streamice_bg_surf_slope_y+sy) * 
                0336 !      &         H_streamice(i,j,bi,bj) * rA(i,j,bi,bj)
                0337 !               
                0338 !              ENDIF
                0339 !             ENDDO
                0340 !            ENDDO
                0341 
                0342            IF (float_frac_streamice(i,j,bi,bj) .eq. 1.0) then
                0343 #ifdef USE_ALT_RLOW
                0344             neu_val = .5 * gravity * 
                0345      &       (streamice_density * H_streamice (i,j,bi,bj) ** 2 - 
                0346      &        streamice_density_ocean_avg * R_low_si(i,j,bi,bj) ** 2)
                0347 #else
                0348             neu_val = .5 * gravity * 
                0349      &       (streamice_density * H_streamice (i,j,bi,bj) ** 2 - 
                0350      &        streamice_density_ocean_avg * R_low(i,j,bi,bj) ** 2)
                0351 #endif
                0352            ELSE
                0353             neu_val = .5 * gravity * 
                0354      &       (1-streamice_density/streamice_density_ocean_avg) * 
                0355      &        streamice_density * H_streamice(i,j,bi,bj) ** 2 
                0356            ENDIF
                0357 
                0358            IF ((STREAMICE_ufacemask(i,j,bi,bj) .eq. 2) 
                0359      &      .OR. (STREAMICE_hmask(i-1,j,bi,bj) .eq. 0)
                0360      &      .OR. (STREAMICE_hmask(i-1,j,bi,bj) .eq. 2) ) THEN ! left face of the cell is at a stress boundary
                0361           ! the depth-integrated longitudinal stress is equal to the difference of depth-integrated pressure on either side of the face
                0362           ! on the ice side, it is rho g h^2 / 2
                0363           ! on the ocean side, it is rhow g (delta OD)^2 / 2
                0364           ! OD can be zero under the ice; but it is ASSUMED on the ice-free side of the face, topography elevation is not above the base of the 
                0365           !     ice in the current cell
                0366              
                0367              taudx_SI(i,j,bi,bj) = taudx_SI(i,j,bi,bj) - 
                0368      &        .5 * dyG(i,j,bi,bj)*(neu_val+streamice_addl_backstress)
                0369              taudx_SI(i,j+1,bi,bj) = taudx_SI(i,j+1,bi,bj) - 
                0370      &        .5 * dyG(i,j,bi,bj)*(neu_val+streamice_addl_backstress)
                0371            ENDIF
                0372 
                0373            IF ((STREAMICE_ufacemask(i+1,j,bi,bj) .eq. 2) 
                0374      &      .OR. (STREAMICE_hmask(i+1,j,bi,bj) .eq. 0)
                0375      &      .OR. (STREAMICE_hmask(i+1,j,bi,bj) .eq. 2) ) THEN 
                0376              
                0377              taudx_SI(i+1,j,bi,bj) = taudx_SI(i+1,j,bi,bj) + 
                0378      &        .5 * dyG(i+1,j,bi,bj)*(neu_val+streamice_addl_backstress)  ! note negative sign is due to direction of normal vector
                0379              taudx_SI(i+1,j+1,bi,bj) = taudx_SI(i+1,j+1,bi,bj) +
                0380      &        .5 * dyG(i+1,j,bi,bj)*(neu_val+streamice_addl_backstress)
                0381            ENDIF
                0382 
                0383            IF ((STREAMICE_vfacemask(i,j,bi,bj) .eq. 2) 
                0384      &      .OR. (STREAMICE_hmask(i,j-1,bi,bj) .eq. 0)
                0385      &      .OR. (STREAMICE_hmask(i,j-1,bi,bj) .eq. 2) ) THEN 
                0386              
                0387              taudy_SI(i,j,bi,bj) = taudy_SI(i,j,bi,bj) - 
                0388      &        .5 * dxG(i,j,bi,bj)*(neu_val+streamice_addl_backstress)
                0389              taudy_SI(i+1,j,bi,bj) = taudy_SI(i+1,j,bi,bj) - 
                0390      &        .5 * dxG(i,j,bi,bj)*(neu_val+streamice_addl_backstress)
                0391            ENDIF
                0392 
                0393            IF ((STREAMICE_vfacemask(i,j+1,bi,bj) .eq. 2) 
                0394      &      .OR. (STREAMICE_hmask(i,j+1,bi,bj) .eq. 0)
                0395      &      .OR. (STREAMICE_hmask(i,j+1,bi,bj) .eq. 2) ) THEN 
                0396              
                0397              taudy_SI(i,j+1,bi,bj) = taudy_SI(i,j+1,bi,bj) +
                0398      &        .5 * dxG(i,j+1,bi,bj)*(neu_val+streamice_addl_backstress)
                0399              taudy_SI(i+1,j+1,bi,bj) = taudy_SI(i+1,j+1,bi,bj) +
                0400      &        .5 * dxG(i,j+1,bi,bj)*(neu_val+streamice_addl_backstress)
                0401            ENDIF
                0402 
                0403           ENDIF
                0404          ENDDO
                0405         ENDDO
                0406        ENDDO
                0407       ENDDO
                0408 
                0409 
                0410 
                0411 #endif
                0412       RETURN
                0413       END
                0414        
                0415
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!