C $Header: /u/gcmpack/MITgcm/pkg/streamice/streamice_adv_front.F,v 1.5 2016/10/24 14:13:12 dgoldberg Exp $
C $Name: $
#include "STREAMICE_OPTIONS.h"
#ifdef ALLOW_AUTODIFF
# include "AUTODIFF_OPTIONS.h"
#endif
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
CBOP
SUBROUTINE STREAMICE_ADV_FRONT (
& myThid,
& time_step,
& hflux_x_si,
& hflux_y_si )
C /============================================================\
C | SUBROUTINE |
C | o |
C |============================================================|
C | |
C \============================================================/
IMPLICIT NONE
C === Global variables ===
#include "SIZE.h"
#include "GRID.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "STREAMICE.h"
#include "STREAMICE_ADV.h"
#ifdef ALLOW_AUTODIFF_TAMC
# include "tamc.h"
#endif
INTEGER myThid
_RL time_step
_RL hflux_x_SI (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
_RL hflux_y_SI (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
#ifdef ALLOW_STREAMICE
INTEGER i, j, bi, bj, k, iter_count, iter_rpt
INTEGER Gi, Gj
INTEGER new_partial(4)
INTEGER ikey_front, ikey_1
_RL iter_flag
_RL n_flux_1, n_flux_2
_RL href, rho, partial_vol, tot_flux, hpot
CHARACTER*(MAX_LEN_MBUF) msgBuf
_RL hflux_x_SI2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
_RL hflux_y_SI2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
rho = streamice_density
cph iter_count = 0
iter_flag = 1. _d 0
iter_rpt = 0
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
hflux_x_SI2(i,j,bi,bj) = 0. _d 0
hflux_y_SI2(i,j,bi,bj) = 0. _d 0
ENDDO
ENDDO
ENDDO
ENDDO
DO iter_count = 0, 3
#ifdef ALLOW_AUTODIFF_TAMC
ikey_front = (ikey_dynamics-1)*4 + iter_count + 1
CADJ STORE area_shelf_streamice
CADJ & = comlev1_stream_front, key = ikey_front
CADJ STORE h_streamice
CADJ & = comlev1_stream_front, key = ikey_front
CADJ STORE hflux_x_si
CADJ & = comlev1_stream_front, key = ikey_front
CADJ STORE hflux_x_si2
CADJ & = comlev1_stream_front, key = ikey_front
CADJ STORE hflux_y_si
CADJ & = comlev1_stream_front, key = ikey_front
CADJ STORE hflux_y_si2
CADJ & = comlev1_stream_front, key = ikey_front
CADJ STORE streamice_hmask
CADJ & = comlev1_stream_front, key = ikey_front
CADJ STORE iter_flag
CADJ & = comlev1_stream_front, key = ikey_front
#endif
IF ( iter_flag .GT. 0. ) THEN
iter_flag = 0. _d 0
IF (iter_count .gt. 0) then
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
hflux_x_SI(i,j,bi,bj)=hflux_x_SI2(i,j,bi,bj)
hflux_y_SI(i,j,bi,bj)=hflux_y_SI2(i,j,bi,bj)
hflux_x_SI2(i,j,bi,bj) = 0. _d 0
hflux_y_SI2(i,j,bi,bj) = 0. _d 0
ENDDO
ENDDO
ENDDO
ENDDO
ENDIF
! iter_count = iter_count + 1
iter_rpt = iter_rpt + 1
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
DO j=1-1,sNy+1
Gj = (myYGlobalLo-1)+(bj-1)*sNy+j
IF ((Gj .ge. 1) .and. (Gj .le. Ny)) THEN
DO i=1-1,sNx+1
Gi = (myXGlobalLo-1)+(bi-1)*sNx+i
#ifdef ALLOW_AUTODIFF_TAMC
act1 = bi - myBxLo(myThid)
max1 = myBxHi(myThid) - myBxLo(myThid) + 1
act2 = bj - myByLo(myThid)
max2 = myByHi(myThid) - myByLo(myThid) + 1
act3 = myThid - 1
max3 = nTx*nTy
act4 = ikey_front - 1
ikey_1 = i + 1
& + (sNx+2)*(j)
& + (sNx+2)*(sNy+2)*act1
& + (sNx+2)*(sNy+2)*max1*act2
& + (sNx+2)*(sNy+2)*max1*max2*act3
& + (sNx+2)*(sNy+2)*max1*max2*max3*act4
CADJ STORE area_shelf_streamice(i,j,bi,bj)
CADJ & = comlev1_stream_ij, key = ikey_1
CADJ STORE h_streamice(i,j,bi,bj)
CADJ & = comlev1_stream_ij, key = ikey_1
CADJ STORE hflux_x_si(i,j,bi,bj)
CADJ & = comlev1_stream_ij, key = ikey_1
CADJ STORE hflux_y_si(i,j,bi,bj)
CADJ & = comlev1_stream_ij, key = ikey_1
CADJ STORE streamice_hmask(i,j,bi,bj)
CADJ & = comlev1_stream_ij, key = ikey_1
#endif
IF (.not. STREAMICE_calve_to_mask .OR.
& STREAMICE_calve_mask (i,j,bi,bj) .eq. 1.0) THEN
IF ((Gi .ge. 1) .and. (Gi .le. Nx) .and.
& (STREAMICE_Hmask(i,j,bi,bj).eq.0.0 .or.
& STREAMICE_Hmask(i,j,bi,bj).eq.2.0)) THEN
n_flux_1 = 0. _d 0
href = 0. _d 0
tot_flux = 0. _d 0
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE hflux_x_SI(i,j,bi,bj)
CADJ & = comlev1_stream_ij, key = ikey_1
#endif
IF (hflux_x_SI(i,j,bi,bj).gt. 0. _d 0) THEN
n_flux_1 = n_flux_1 + 1. _d 0
href = href + H_streamice(i-1,j,bi,bj)
tot_flux = tot_flux + hflux_x_SI(i,j,bi,bj) *
& dxG(i,j,bi,bj) * time_step
hflux_x_SI(i,j,bi,bj) = 0. _d 0
ENDIF
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE hflux_x_SI(i,j,bi,bj)
CADJ & = comlev1_stream_ij, key = ikey_1
#endif
IF (hflux_x_SI(i+1,j,bi,bj).lt. 0. _d 0) THEN
n_flux_1 = n_flux_1 + 1. _d 0
href = href + H_streamice(i+1,j,bi,bj)
tot_flux = tot_flux - hflux_x_SI(i+1,j,bi,bj) *
& dxG(i+1,j,bi,bj) * time_step
hflux_x_SI(i+1,j,bi,bj) = 0. _d 0
ENDIF
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE hflux_y_SI(i,j,bi,bj)
CADJ & = comlev1_stream_ij, key = ikey_1
#endif
IF (hflux_y_SI(i,j,bi,bj).gt. 0. _d 0) THEN
n_flux_1 = n_flux_1 + 1. _d 0
href = href + H_streamice(i,j-1,bi,bj)
tot_flux = tot_flux + hflux_y_SI(i,j,bi,bj) *
& dyG(i,j,bi,bj) * time_step
hflux_y_SI(i,j,bi,bj) = 0. _d 0
ENDIF
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE hflux_y_SI(i,j,bi,bj)
CADJ & = comlev1_stream_ij, key = ikey_1
#endif
IF (hflux_y_SI(i,j+1,bi,bj).lt. 0. _d 0) THEN
n_flux_1 = n_flux_1 + 1. _d 0
href = href + H_streamice(i,j+1,bi,bj)
tot_flux = tot_flux - hflux_y_SI(i,j+1,bi,bj) *
& dyG(i,j+1,bi,bj) * time_step
hflux_y_SI(i,j+1,bi,bj) = 0. _d 0
ENDIF
IF (n_flux_1 .gt. 0.) THEN
href = href / n_flux_1
partial_vol = H_streamice (i,j,bi,bj) *
& area_shelf_streamice (i,j,bi,bj) + tot_flux
hpot = partial_vol * recip_rA(i,j,bi,bj)
IF (hpot .eq. href) THEN ! cell is exactly covered, no overflow
STREAMICE_hmask (i,j,bi,bj) = 1.0
H_streamice (i,j,bi,bj) = href
area_shelf_streamice(i,j,bi,bj) =
& rA(i,j,bi,bj)
ELSEIF (hpot .lt. href) THEN ! cell still unfilled
STREAMICE_hmask (i,j,bi,bj) = 2.0
area_shelf_streamice (i,j,bi,bj) = partial_vol / href
H_streamice (i,j,bi,bj) = href
ELSE ! cell is filled - do overflow
STREAMICE_hmask (i,j,bi,bj) = 1.0
area_shelf_streamice(i,j,bi,bj) =
& rA(i,j,bi,bj)
PRINT *, "GOT HERE OVERFLOW ", i,j,
& area_shelf_streamice(i,j,bi,bj)
partial_vol = partial_vol - href * rA(i,j,bi,bj)
iter_flag = 1. _d 0
n_flux_2 = 0. _d 0
DO k=1,4
new_partial (:) = 0
ENDDO
DO k=1,2
IF ( (STREAMICE_ufacemask(i-1+k,j,bi,bj).eq.2.0) .or.
& (STREAMICE_calve_to_mask .and.
& STREAMICE_calve_mask(i+2*k-3,j,bi,bj).ne.1.0)
& ) THEN ! at a permanent calving boundary - no advance allowed
n_flux_2 = n_flux_2 + 1. _d 0
ELSEIF (STREAMICE_hmask(i+2*k-3,j,bi,bj).eq.0 _d 0) THEN ! adjacent cell is completely ice free
n_flux_2 = n_flux_2 + 1. _d 0
new_partial (k) = 1
ENDIF
ENDDO
DO k=1,2
IF ( (STREAMICE_vfacemask (i,j-1+k,bi,bj).eq.2.0) .or.
& (STREAMICE_calve_to_mask .and.
& STREAMICE_calve_mask(i,j+2*k-3,bi,bj).ne.1.0)
& ) THEN ! at a permanent calving boundary - no advance allowed
n_flux_2 = n_flux_2 + 1. _d 0
ELSEIF (STREAMICE_hmask(i,j+2*k-3,bi,bj).eq.0 _d 0) THEN
n_flux_2 = n_flux_2 + 1. _d 0
new_partial (k+2) = 1
ENDIF
ENDDO
IF (n_flux_2 .eq. 0.) THEN ! there is nowhere to put the extra ice!
H_streamice(i,j,bi,bj) = href + partial_vol *
& recip_rA(i,j,bi,bj)
ELSE
H_streamice(i,j,bi,bj) = href
DO k=1,2
IF (new_partial(k) .eq. 1) THEN
hflux_x_SI2(i-1+k,j,bi,bj) =
& partial_vol/time_step/n_flux_2/
& dxG(i-1+k,j,bi,bj)
ENDIF
ENDDO
DO k=1,2
IF (new_partial(k+2) .eq. 1) THEN
hflux_y_SI2(i,j-1+k,bi,bj) =
& partial_vol/time_step/n_flux_2/
& dxG(i,j-1+k,bi,bj)
ENDIF
ENDDO
ENDIF
ENDIF
ENDIF
ENDIF
ENDIF
ENDDO
ENDIF
ENDDO
c
ENDDO
ENDDO
c
ENDIF
ENDDO
IF (iter_rpt.gt.1) THEN
WRITE(msgBuf,'(A,I5,A)') 'FRONT ADVANCE: ',iter_rpt,
& ' ITERATIONS'
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
& SQUEEZE_RIGHT , 1)
ENDIF
#endif
RETURN
END