Back to home page

MITgcm

 
 

    


File indexing completed on 2023-09-21 05:10:47 UTC

view on githubraw file Latest commit 96b00645 on 2023-09-20 15:15:14 UTC
5ca83cd8f7 Dani*0001 #include "STREAMICE_OPTIONS.h"
9d8b0af494 Jean*0002 #ifdef ALLOW_AUTODIFF
                0003 # include "AUTODIFF_OPTIONS.h"
                0004 #endif
5ca83cd8f7 Dani*0005 
                0006 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0007 
                0008 CBOP
9d8b0af494 Jean*0009       SUBROUTINE STREAMICE_ADV_FRONT (
                0010      & myThid,
5ca83cd8f7 Dani*0011      & time_step,
                0012      & hflux_x_si,
                0013      & hflux_y_si )
                0014 
                0015 C     /============================================================\
9d8b0af494 Jean*0016 C     | SUBROUTINE                                                 |
5ca83cd8f7 Dani*0017 C     | o                                                          |
                0018 C     |============================================================|
                0019 C     |                                                            |
                0020 C     \============================================================/
                0021       IMPLICIT NONE
                0022 
                0023 C     === Global variables ===
                0024 #include "SIZE.h"
                0025 #include "GRID.h"
                0026 #include "EEPARAMS.h"
                0027 #include "PARAMS.h"
                0028 #include "STREAMICE.h"
                0029 #include "STREAMICE_ADV.h"
                0030 #ifdef ALLOW_AUTODIFF_TAMC
                0031 # include "tamc.h"
                0032 #endif
                0033 
                0034       INTEGER myThid
                0035       _RL time_step
                0036       _RL hflux_x_SI (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0037       _RL hflux_y_SI (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0038 
                0039 #ifdef ALLOW_STREAMICE
                0040 
                0041       INTEGER i, j, bi, bj, k, iter_count, iter_rpt
                0042       INTEGER Gi, Gj
                0043       INTEGER new_partial(4)
7c50f07931 Mart*0044 #ifdef ALLOW_AUTODIFF_TAMC
5ca83cd8f7 Dani*0045       INTEGER ikey_front, ikey_1
7c50f07931 Mart*0046 #endif
5ca83cd8f7 Dani*0047       _RL iter_flag
                0048       _RL n_flux_1, n_flux_2
                0049       _RL href, rho, partial_vol, tot_flux, hpot
                0050       CHARACTER*(MAX_LEN_MBUF) msgBuf
                0051       _RL hflux_x_SI2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0052       _RL hflux_y_SI2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0053 
                0054       rho = streamice_density
                0055 cph      iter_count = 0
                0056       iter_flag = 1. _d 0
                0057       iter_rpt = 0
                0058 
                0059         DO bj=myByLo(myThid),myByHi(myThid)
                0060          DO bi=myBxLo(myThid),myBxHi(myThid)
                0061           DO j=1-OLy,sNy+OLy
                0062            DO i=1-OLx,sNx+OLx
                0063             hflux_x_SI2(i,j,bi,bj) = 0. _d 0
                0064             hflux_y_SI2(i,j,bi,bj) = 0. _d 0
                0065            ENDDO
                0066           ENDDO
                0067          ENDDO
                0068         ENDDO
                0069 
                0070       DO iter_count = 0, 3
                0071 
                0072 #ifdef ALLOW_AUTODIFF_TAMC
                0073          ikey_front = (ikey_dynamics-1)*4 + iter_count + 1
9d8b0af494 Jean*0074 CADJ STORE area_shelf_streamice
5ca83cd8f7 Dani*0075 CADJ &     = comlev1_stream_front, key = ikey_front
                0076 CADJ STORE h_streamice
                0077 CADJ &     = comlev1_stream_front, key = ikey_front
                0078 CADJ STORE hflux_x_si
                0079 CADJ &     = comlev1_stream_front, key = ikey_front
                0080 CADJ STORE hflux_x_si2
                0081 CADJ &     = comlev1_stream_front, key = ikey_front
                0082 CADJ STORE hflux_y_si
                0083 CADJ &     = comlev1_stream_front, key = ikey_front
                0084 CADJ STORE hflux_y_si2
                0085 CADJ &     = comlev1_stream_front, key = ikey_front
                0086 CADJ STORE streamice_hmask
                0087 CADJ &     = comlev1_stream_front, key = ikey_front
                0088 CADJ STORE iter_flag
                0089 CADJ &     = comlev1_stream_front, key = ikey_front
                0090 #endif
                0091 
                0092        IF ( iter_flag .GT. 0. ) THEN
                0093 
                0094        iter_flag = 0. _d 0
                0095 
                0096        IF (iter_count .gt. 0) then
                0097         DO bj=myByLo(myThid),myByHi(myThid)
                0098          DO bi=myBxLo(myThid),myBxHi(myThid)
                0099           DO j=1-OLy,sNy+OLy
                0100            DO i=1-OLx,sNx+OLx
                0101             hflux_x_SI(i,j,bi,bj)=hflux_x_SI2(i,j,bi,bj)
                0102             hflux_y_SI(i,j,bi,bj)=hflux_y_SI2(i,j,bi,bj)
                0103             hflux_x_SI2(i,j,bi,bj) = 0. _d 0
                0104             hflux_y_SI2(i,j,bi,bj) = 0. _d 0
                0105            ENDDO
                0106           ENDDO
                0107          ENDDO
9d8b0af494 Jean*0108         ENDDO
5ca83cd8f7 Dani*0109        ENDIF
                0110 
96b006450c dngo*0111 c       iter_count = iter_count + 1
5ca83cd8f7 Dani*0112        iter_rpt = iter_rpt + 1
                0113 
                0114        DO bj=myByLo(myThid),myByHi(myThid)
                0115         DO bi=myBxLo(myThid),myBxHi(myThid)
                0116 
                0117          DO j=1-1,sNy+1
                0118           Gj = (myYGlobalLo-1)+(bj-1)*sNy+j
                0119           IF ((Gj .ge. 1) .and. (Gj .le. Ny)) THEN
                0120            DO i=1-1,sNx+1
                0121             Gi = (myXGlobalLo-1)+(bi-1)*sNx+i
                0122 
                0123 #ifdef ALLOW_AUTODIFF_TAMC
20dee61641 Mart*0124             ikey_1 = bi + (bj-1)*nSx + (ikey_front-1)*nSx*nSy
                0125             ikey_1 = i + 1  + (j)*(sNx+2) + (ikey_1-1)*(sNx+2)*(sNy+2)
5ca83cd8f7 Dani*0126 CADJ STORE area_shelf_streamice(i,j,bi,bj)
                0127 CADJ &     = comlev1_stream_ij, key = ikey_1
                0128 CADJ STORE h_streamice(i,j,bi,bj)
                0129 CADJ &     = comlev1_stream_ij, key = ikey_1
                0130 CADJ STORE hflux_x_si(i,j,bi,bj)
                0131 CADJ &     = comlev1_stream_ij, key = ikey_1
                0132 CADJ STORE hflux_y_si(i,j,bi,bj)
                0133 CADJ &     = comlev1_stream_ij, key = ikey_1
                0134 CADJ STORE streamice_hmask(i,j,bi,bj)
                0135 CADJ &     = comlev1_stream_ij, key = ikey_1
                0136 #endif
                0137 
a84d06a708 Dani*0138             IF (.not. STREAMICE_calve_to_mask .OR.
eaf63fbcc2 Dani*0139      &       STREAMICE_calve_mask (i,j,bi,bj) .eq. 1.0) THEN
                0140 
5ca83cd8f7 Dani*0141             IF ((Gi .ge. 1) .and. (Gi .le. Nx) .and.
                0142      &          (STREAMICE_Hmask(i,j,bi,bj).eq.0.0 .or.
                0143      &           STREAMICE_Hmask(i,j,bi,bj).eq.2.0)) THEN
                0144              n_flux_1 = 0. _d 0
                0145              href = 0. _d 0
                0146              tot_flux = 0. _d 0
                0147 
                0148 #ifdef ALLOW_AUTODIFF_TAMC
96b006450c dngo*0149 CADJ STORE hflux_x_SI(i,j,bi,bj) = comlev1_stream_ij, key = ikey_1
5ca83cd8f7 Dani*0150 #endif
                0151              IF (hflux_x_SI(i,j,bi,bj).gt. 0. _d 0) THEN
                0152               n_flux_1 = n_flux_1 + 1. _d 0
                0153               href = href + H_streamice(i-1,j,bi,bj)
9d8b0af494 Jean*0154               tot_flux = tot_flux + hflux_x_SI(i,j,bi,bj) *
5ca83cd8f7 Dani*0155      &         dxG(i,j,bi,bj) * time_step
                0156               hflux_x_SI(i,j,bi,bj) = 0. _d 0
                0157              ENDIF
                0158 
                0159 #ifdef ALLOW_AUTODIFF_TAMC
96b006450c dngo*0160 CADJ STORE hflux_x_SI(i+1,j,bi,bj) = comlev1_stream_ij, key = ikey_1
5ca83cd8f7 Dani*0161 #endif
                0162              IF (hflux_x_SI(i+1,j,bi,bj).lt. 0. _d 0) THEN
                0163               n_flux_1 = n_flux_1 + 1. _d 0
                0164               href = href + H_streamice(i+1,j,bi,bj)
9d8b0af494 Jean*0165               tot_flux = tot_flux - hflux_x_SI(i+1,j,bi,bj) *
5ca83cd8f7 Dani*0166      &         dxG(i+1,j,bi,bj) * time_step
                0167               hflux_x_SI(i+1,j,bi,bj) = 0. _d 0
                0168              ENDIF
                0169 
                0170 #ifdef ALLOW_AUTODIFF_TAMC
96b006450c dngo*0171 CADJ STORE hflux_y_SI(i,j,bi,bj) = comlev1_stream_ij, key = ikey_1
5ca83cd8f7 Dani*0172 #endif
                0173              IF (hflux_y_SI(i,j,bi,bj).gt. 0. _d 0) THEN
                0174               n_flux_1 = n_flux_1 + 1. _d 0
                0175               href = href + H_streamice(i,j-1,bi,bj)
9d8b0af494 Jean*0176               tot_flux = tot_flux + hflux_y_SI(i,j,bi,bj) *
5ca83cd8f7 Dani*0177      &         dyG(i,j,bi,bj) * time_step
                0178               hflux_y_SI(i,j,bi,bj) = 0. _d 0
                0179              ENDIF
                0180 
                0181 #ifdef ALLOW_AUTODIFF_TAMC
96b006450c dngo*0182 CADJ STORE hflux_y_SI(i,j+1,bi,bj) = comlev1_stream_ij, key = ikey_1
5ca83cd8f7 Dani*0183 #endif
                0184              IF (hflux_y_SI(i,j+1,bi,bj).lt. 0. _d 0) THEN
                0185               n_flux_1 = n_flux_1 + 1. _d 0
                0186               href = href + H_streamice(i,j+1,bi,bj)
9d8b0af494 Jean*0187               tot_flux = tot_flux - hflux_y_SI(i,j+1,bi,bj) *
5ca83cd8f7 Dani*0188      &         dyG(i,j+1,bi,bj) * time_step
                0189               hflux_y_SI(i,j+1,bi,bj) = 0. _d 0
                0190              ENDIF
9d8b0af494 Jean*0191 
5ca83cd8f7 Dani*0192              IF (n_flux_1 .gt. 0.) THEN
                0193 
                0194               href = href / n_flux_1
9d8b0af494 Jean*0195               partial_vol = H_streamice (i,j,bi,bj) *
5ca83cd8f7 Dani*0196      &         area_shelf_streamice (i,j,bi,bj) + tot_flux
                0197               hpot = partial_vol * recip_rA(i,j,bi,bj)
9d8b0af494 Jean*0198 
5ca83cd8f7 Dani*0199               IF (hpot .eq. href) THEN ! cell is exactly covered, no overflow
                0200                 STREAMICE_hmask (i,j,bi,bj) = 1.0
                0201                 H_streamice (i,j,bi,bj) = href
9d8b0af494 Jean*0202                 area_shelf_streamice(i,j,bi,bj) =
5ca83cd8f7 Dani*0203      &           rA(i,j,bi,bj)
                0204               ELSEIF (hpot .lt. href) THEN ! cell still unfilled
                0205 
                0206                STREAMICE_hmask (i,j,bi,bj) = 2.0
                0207                area_shelf_streamice (i,j,bi,bj) = partial_vol / href
                0208                H_streamice (i,j,bi,bj) = href
                0209               ELSE ! cell is filled - do overflow
                0210 
                0211                STREAMICE_hmask (i,j,bi,bj) = 1.0
9d8b0af494 Jean*0212                area_shelf_streamice(i,j,bi,bj) =
5ca83cd8f7 Dani*0213      &           rA(i,j,bi,bj)
                0214 
9d8b0af494 Jean*0215                PRINT *, "GOT HERE OVERFLOW ", i,j,
eaf63fbcc2 Dani*0216      &          area_shelf_streamice(i,j,bi,bj)
5ca83cd8f7 Dani*0217                partial_vol = partial_vol - href * rA(i,j,bi,bj)
                0218 
                0219                iter_flag  = 1. _d 0
                0220 
9d8b0af494 Jean*0221                n_flux_2 = 0. _d 0 ;
5ca83cd8f7 Dani*0222                DO k=1,4
                0223                 new_partial (:) = 0
                0224                ENDDO
                0225 
                0226                DO k=1,2
9d8b0af494 Jean*0227                 IF ( (STREAMICE_ufacemask(i-1+k,j,bi,bj).eq.2.0) .or.
                0228      &            (STREAMICE_calve_to_mask .and.
                0229      &             STREAMICE_calve_mask(i+2*k-3,j,bi,bj).ne.1.0)
eaf63fbcc2 Dani*0230      &             ) THEN  ! at a permanent calving boundary - no advance allowed
5ca83cd8f7 Dani*0231                    n_flux_2 = n_flux_2 + 1. _d 0
                0232                 ELSEIF (STREAMICE_hmask(i+2*k-3,j,bi,bj).eq.0 _d 0) THEN ! adjacent cell is completely ice free
                0233                    n_flux_2 = n_flux_2 + 1. _d 0
                0234                    new_partial (k) = 1
                0235                 ENDIF
                0236                ENDDO
                0237                DO k=1,2
9d8b0af494 Jean*0238                 IF ( (STREAMICE_vfacemask (i,j-1+k,bi,bj).eq.2.0) .or.
                0239      &            (STREAMICE_calve_to_mask .and.
                0240      &             STREAMICE_calve_mask(i,j+2*k-3,bi,bj).ne.1.0)
eaf63fbcc2 Dani*0241      &             ) THEN  ! at a permanent calving boundary - no advance allowed
5ca83cd8f7 Dani*0242                     n_flux_2 = n_flux_2 + 1. _d 0
                0243                 ELSEIF (STREAMICE_hmask(i,j+2*k-3,bi,bj).eq.0 _d 0) THEN
                0244                     n_flux_2 = n_flux_2 + 1. _d 0
                0245                     new_partial (k+2) = 1
                0246                 ENDIF
                0247                ENDDO
                0248 
                0249                IF (n_flux_2 .eq. 0.) THEN ! there is nowhere to put the extra ice!
                0250                 H_streamice(i,j,bi,bj) = href + partial_vol *
                0251      &             recip_rA(i,j,bi,bj)
                0252                ELSE
                0253                 H_streamice(i,j,bi,bj) = href
                0254 
                0255                 DO k=1,2
                0256                  IF (new_partial(k) .eq. 1) THEN
9d8b0af494 Jean*0257                   hflux_x_SI2(i-1+k,j,bi,bj) =
5ca83cd8f7 Dani*0258      &             partial_vol/time_step/n_flux_2/
                0259      &               dxG(i-1+k,j,bi,bj)
                0260                  ENDIF
                0261                 ENDDO
                0262 
                0263                 DO k=1,2
                0264                  IF (new_partial(k+2) .eq. 1) THEN
9d8b0af494 Jean*0265                   hflux_y_SI2(i,j-1+k,bi,bj) =
5ca83cd8f7 Dani*0266      &             partial_vol/time_step/n_flux_2/
                0267      &               dxG(i,j-1+k,bi,bj)
                0268                  ENDIF
                0269                 ENDDO
                0270 
                0271                ENDIF
                0272               ENDIF
                0273              ENDIF
                0274 
                0275             ENDIF
eaf63fbcc2 Dani*0276             ENDIF
5ca83cd8f7 Dani*0277            ENDDO
                0278           ENDIF
                0279          ENDDO
                0280 c
                0281         ENDDO
                0282        ENDDO
                0283 c
                0284       ENDIF
                0285       ENDDO
                0286 
                0287       IF (iter_rpt.gt.1) THEN
                0288        WRITE(msgBuf,'(A,I5,A)') 'FRONT ADVANCE: ',iter_rpt,
                0289      &  ' ITERATIONS'
                0290        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
9d8b0af494 Jean*0291      &                     SQUEEZE_RIGHT , 1)
5ca83cd8f7 Dani*0292       ENDIF
                0293 
                0294 #endif
                0295       RETURN
                0296       END