MODULE ICEDYN_COM 9,2
!@sum ICEDYN_COM holds global variables for dynamic sea ice
!@auth Gary Russell/Gavin Schmidt
!@ver 1.0
USE MODEL_COM
, only : im,jm
USE ICEDYN, only : imic,jmic
IMPLICIT NONE
SAVE
C**** variables used in advection (on ICE grid)
!@var RSIX,RSIY first order moments for seaice concentration
!@var USI,VSI east-west, and north-south sea ice velocities (m/s)
! REAL*8, DIMENSION(IMIC,JMIC) :: RSIX,RSIY,USI,VSI
REAL*8, ALLOCATABLE, DIMENSION(:,:) :: RSIX,RSIY,USI,VSI
!@var USIDT,VSIDT sea ice fluxes, saved for advection (m)
! REAL*8, DIMENSION(IMIC,JMIC) :: USIDT,VSIDT
REAL*8, ALLOCATABLE, DIMENSION(:,:) :: USIDT,VSIDT
C**** Needed for ADVSI (on ATM grid)
!@var RSISAVE saved value of sea ice concentration before DYNSI
! REAL*8, DIMENSION(IM,JM) :: RSISAVE
REAL*8, ALLOCATABLE, DIMENSION(:,:) :: RSISAVE
C**** Ice advection diagnostics
INTEGER, PARAMETER :: KICIJ=12
!@var IJ_xxx Names for ICIJ diagnostics
INTEGER IJ_USI,IJ_VSI,IJ_DMUI,IJ_DMVI,IJ_PICE,IJ_MUSI,IJ_MVSI
* ,IJ_HUSI,IJ_HVSI,IJ_SUSI,IJ_SVSI,IJ_RSI
!@var ICIJ lat-lon ice dynamic diagnostics (on atm grid)
! REAL*8, DIMENSION(IMIC,JMIC,KICIJ) :: ICIJ
REAL*8, ALLOCATABLE, DIMENSION(:,:,:) :: ICIJ
!@var lname_icij Long names for ICIJ diagnostics
CHARACTER*50, DIMENSION(KICIJ) :: LNAME_ICIJ
!@var sname_icij Short names for ICIJ diagnostics
CHARACTER*30, DIMENSION(KICIJ) :: SNAME_ICIJ
!@var units_icij Units for ICIJ diagnostics
CHARACTER*50, DIMENSION(KICIJ) :: UNITS_ICIJ
!@var ia_icij IDACC numbers for ICIJ diagnostics
INTEGER, DIMENSION(KICIJ) :: IA_ICIJ
!@var scale_icij scales for ICIJ diagnostics
REAL*8, DIMENSION(KICIJ) :: SCALE_ICIJ
!@var ijgrid_icij Grid descriptor for ICIJ diagnostics
INTEGER, DIMENSION(KICIJ) :: IJGRID_ICIJ
END MODULE ICEDYN_COM
SUBROUTINE ALLOC_ICEDYN_COM(grid) 1,9
!@sum ALLOC_ICEDYN_COM allocates arrays defined in the ICEDYN_COM module
!@auth Rosalinda de Fainchtein
USE DOMAIN_DECOMP, only : GET
USE DOMAIN_DECOMP, only : DYN_GRID
USE MODEL_COM, only : im
USE ICEDYN, only : imic
USE ICEDYN, only : grid_MIC
USE ICEDYN_COM, only : KICIJ
USE ICEDYN_COM, only : RSIX,RSIY,USI,VSI,USIDT,VSIDT,
& RSISAVE,ICIJ
IMPLICIT NONE
LOGICAL, SAVE :: init=.false.
INTEGER :: J_1H , J_0H
INTEGER :: J_1H_MIC, J_0H_MIC
INTEGER :: IER
TYPE(DYN_GRID) :: grid
If (init) Then
Return ! Only invoke once
End If
init = .true.
!*** For now set grid_MIC to be the same as grid
! This is consistent with the current status of the code for
! parallelization along latitude (j)
grid_MIC =grid
C**** Get dimensioning parameters for arrays defined in the grid and grid_MIC
C**** stencils.
CALL GET(grid , J_STRT_HALO=J_0H , J_STOP_HALO=J_1H )
CALL GET(grid_MIC, J_STRT_HALO=J_0H_MIC, J_STOP_HALO=J_1H_MIC)
ALLOCATE( RSIX(IMIC, J_0H_MIC:J_1H_MIC),
& RSIY(IMIC, J_0H_MIC:J_1H_MIC),
& USI(IMIC, J_0H_MIC:J_1H_MIC),
& VSI(IMIC, J_0H_MIC:J_1H_MIC),
& STAT = IER)
ALLOCATE( USIDT(IMIC, J_0H_MIC:J_1H_MIC),
& VSIDT(IMIC, J_0H_MIC:J_1H_MIC),
& STAT = IER)
ALLOCATE( RSISAVE(IM, J_0H:J_1H),
& STAT = IER)
ALLOCATE( ICIJ(IMIC, J_0H_MIC:J_1H_MIC, KICIJ),
& STAT = IER)
return
END SUBROUTINE ALLOC_ICEDYN_COM
SUBROUTINE io_icedyn(kunit,iaction,ioerr) 1,2
!@sum io_icedyn reads and writes dynamic ice arrays to file
!@auth Gavin Schmidt
!@ver 1.0
USE MODEL_COM, only : ioread,iowrite,irsfic,irsficno,irsficnt
* ,irerun,lhead
USE ICEDYN_COM
IMPLICIT NONE
INTEGER kunit !@var kunit unit number of read/write
INTEGER iaction !@var iaction flag for reading or writing to file
!@var IOERR 1 (or -1) if there is (or is not) an error in i/o
INTEGER, INTENT(INOUT) :: IOERR
!@var HEADER Character string label for individual records
CHARACTER*80 :: HEADER, MODULE_HEADER = "ICEDYN01"
write(MODULE_HEADER(lhead+1:80),'(a7,i3,a1,i3,a)')
* 'R8 dim(',imic,',',jmic,'):RSIX,RSIY,USI,VSI'
SELECT CASE (IACTION)
CASE (:IOWRITE) ! output to standard restart file
WRITE (kunit,err=10) MODULE_HEADER,RSIX,RSIY,USI,VSI
CASE (IOREAD:) ! input from restart file
SELECT CASE (IACTION)
CASE (IRSFICNO) ! initial conditions (no ocean)
CASE (ioread,irerun,irsfic,irsficnt) ! restarts
READ (kunit,err=10) HEADER,RSIX,RSIY,USI,VSI
IF (HEADER(1:LHEAD).NE.MODULE_HEADER(1:LHEAD)) THEN
PRINT*,"Discrepancy in module version ",HEADER,MODULE_HEADER
GO TO 10
END IF
END SELECT
END SELECT
RETURN
10 IOERR=1
RETURN
C****
END SUBROUTINE io_icedyn
SUBROUTINE io_icdiag(kunit,it,iaction,ioerr) 1,2
!@sum io_icdiag reads and writes ice dynamic diagnostic arrays to file
!@auth Gavin Schmidt
!@ver 1.0
USE MODEL_COM, only : ioread,iowrite,iowrite_mon,iowrite_single
* ,irsfic,irsficnt,irerun,ioread_single,lhead
USE ICEDYN_COM
IMPLICIT NONE
INTEGER kunit !@var kunit unit number of read/write
INTEGER iaction !@var iaction flag for reading or writing to file
!@var IOERR 1 (or -1) if there is (or is not) an error in i/o
INTEGER, INTENT(INOUT) :: IOERR
!@var HEADER Character string label for individual records
CHARACTER*80 :: HEADER, MODULE_HEADER = "ICDIAG01"
!@var it input/ouput value of hour
INTEGER, INTENT(INOUT) :: it
!@var ICIJ4 dummy arrays for reading diag. files
REAL*4, DIMENSION(IMIC,JMIC,KICIJ) :: ICIJ4
write(MODULE_HEADER(lhead+1:80),'(a8,i3,a1,i3,a1,i2,a4)')
* 'R8 ICij(',imic,',',jmic,',',kicij,'),it'
SELECT CASE (IACTION)
CASE (IOWRITE,IOWRITE_MON) ! output to standard restart file
WRITE (kunit,err=10) MODULE_HEADER,ICIJ,it
CASE (IOWRITE_SINGLE) ! output to acc file
MODULE_HEADER(LHEAD+1:LHEAD+2) = 'R4'
WRITE (kunit,err=10) MODULE_HEADER,REAL(ICIJ,KIND=4),it
CASE (IOREAD:) ! input from restart file
SELECT CASE (IACTION)
CASE (ioread_single) ! accumulate diagnostic files
READ (kunit,err=10) HEADER,ICIJ4,it
C**** accumulate diagnostics
ICIJ=ICIJ+ICIJ4
IF (HEADER(1:LHEAD).NE.MODULE_HEADER(1:LHEAD)) THEN
PRINT*,"Discrepancy in module version ",HEADER
* ,MODULE_HEADER
GO TO 10
END IF
CASE (ioread,irerun) ! restarts
READ (kunit,err=10) HEADER,ICIJ,it
IF (HEADER(1:LHEAD).NE.MODULE_HEADER(1:LHEAD)) THEN
PRINT*,"Discrepancy in module version ",HEADER
* ,MODULE_HEADER
GO TO 10
END IF
CASE (IRSFIC) ! initial conditions
READ (kunit)
CASE (IRSFICNT) ! initial conditions (with no tracers)
READ (kunit)
END SELECT
END SELECT
RETURN
10 IOERR=1
RETURN
C****
END SUBROUTINE io_icdiag
SUBROUTINE reset_icdiag 1,1
!@sum reset_icdiag resets ice dynamic diagnostic arrays
!@auth Gavin Schmidt
USE ICEDYN_COM
IMPLICIT NONE
ICIJ=0.
RETURN
END SUBROUTINE reset_icdiag
SUBROUTINE DYNSI 1,15
!@sum DYNSI calculate ice velocites
!@+ Note that the ice velocities are calculated on the ice grid
!@auth Jiping Liu/Gavin Schmidt (based on code from J. Zhang)
!@ver 1.0
USE CONSTANT, only : rhoi,grav,omega,rhows
USE MODEL_COM, only : im,jm,p,ptop,dts=>dtsrc,focean
USE DOMAIN_DECOMP, only : grid, DYN_GRID, GET
USE DOMAIN_DECOMP, only : CHECKSUM, HALO_UPDATE, NORTH, SOUTH
USE GEOM, only : dxyn,dxys,dxyv,dxyp,bydxyp,dxp,dyv,imaxj
USE ICEDYN, only : imic,jmic,nx1,ny1,press,heffm,uvm,dwatn,cor
* ,sinwat,coswat,bydts,sinen,uice,vice,heff,area,gairx,gairy
* ,gwatx,gwaty,pgfub,pgfvb,amass,uicec,vicec,uib,vib,dmu,dmv
USE ICEDYN_COM, only : usi,vsi,usidt,vsidt,rsisave,icij,ij_usi
* ,ij_vsi,ij_dmui,ij_dmvi,ij_pice,ij_rsi
USE FLUXES, only : dmua,dmva,dmui,dmvi,UI2rho,ogeoza,uosurf,vosurf
* ,apress,uisurf,visurf
USE SEAICE, only : ace1i
USE SEAICE_COM, only : rsi,msi,snowi
IMPLICIT NONE
SAVE
C**** intermediate calculation for pressure gradient terms
REAL*8, DIMENSION(IM, grid%J_STRT_HALO:grid%J_STOP_HALO) ::
& PGFU,PGFV
C****
REAL*8, PARAMETER :: BYRHOI=1D0/RHOI
REAL*8 :: hemi
INTEGER I,J,n,k,ip1,im1,l
REAL*8 USINP,DMUINP,duA,dvA
C**** Declare a new grid data type grid_NXY (to handle do j=1,ny1 type
C Loops. Set grid_NXY=grid, consistent with curr. state of the code.
TYPE(DYN_GRID) :: grid_NXY
INTEGER :: J_1NXY, J_0NXY
INTEGER :: J_1NXYS
INTEGER :: J_1 , J_0
INTEGER :: J_1S , J_0S
INTEGER :: J_1STG,J_0STG
grid_NXY=grid
C**** Get loop indices corresponding to grid and grid_NXY structures
CALL GET(grid_NXY, J_STRT=J_0NXY , J_STOP=J_1NXY
& , J_STOP_SKP=J_1NXYS)
call GET(grid , J_STRT=J_0 , J_STOP=J_1 )
call GET(grid , J_STRT_SKP=J_0S , J_STOP_SKP=J_1S)
call GET(grid , J_STRT_STGR=J_0STG, J_STOP_STGR=J_1STG)
C**** Start main loop
C**** Replicate polar boxes
if (grid%HAVE_NORTH_POLE) then
RSI(2:IM,JM)=RSI(1,JM)
MSI(2:IM,JM)=MSI(1,JM)
DMUA(2:IM,JM,2) = DMUA(1,JM,2)
DMVA(2:IM,JM,2) = DMVA(1,JM,2)
end if
C**** save current value of sea ice concentration for ADVSI
C**** RSISAVE is on atmospheric grid
RSISAVE(:,:)=RSI(:,:)
C**** Pressure anomaly at surface APRESS: calculated by sea ice routines
C**** APRESS is on atmospheric grid. We are now no longer using this as
C**** a forcing in the sea ice dynamics because it is partially
C**** included already in the internal ice pressure gradient. The
C**** atmospheric pressure gradient term does not in general produce a
C**** horizontal force in a solid (such as ice).
C**** calculate sea surface tilt on atmospheric C grid
C**** (using OGEOZA on atmospheric grid plus displacement of free
C**** surface due to presence of ice). This is ignored in favour of
C**** geostrophy if osurf_tilt=0.
C**** PGF is an accelaration
if (grid%HAVE_NORTH_POLE) PGFU(1:IM,JM)=0
if (grid%HAVE_SOUTH_POLE) PGFU(1:IM, 1)=0 !RKF
DO J=J_0S,J_1S
I=IM
DO IP1=1,IM
IF(FOCEAN(I,J).gt.0 .and. FOCEAN(IP1,J).gt.0. .and.
* RSI(I,J)+RSI(IP1,J).gt.0.) THEN
c PGFU(I,J)=-((APRESS(IP1,J)-APRESS(I,J))*BYRHOI
c * +OGEOZA(IP1,J)-OGEOZA(I,J))/DXP(J)
PGFU(I,J)=-(OGEOZA(IP1,J)-OGEOZA(I,J)+
* (RSI(IP1,J)*(MSI(IP1,J)+SNOWI(IP1,J)+ACE1I)
* -RSI(I,J)*(MSI(I,J)+SNOWI(I,J)+ACE1I))*GRAV/RHOWS )
* /DXP(J)
ELSE
PGFU(I,J)=0.
END IF
I=IP1
END DO
END DO
C**** Fill halos for arrays FOCEAN, RSI,OGEOZA,DYV,MSI,SNOWI
C**** Commented Halo fill for array APRESS supports commented statement
CALL CHECKSUM( grid, FOCEAN, 379, "ICEDYN_DRV.f")
CALL HALO_UPDATE(grid, FOCEAN, from=NORTH )
CALL CHECKSUM( grid, RSI , 381, "ICEDYN_DRV.f")
CALL HALO_UPDATE(grid, RSI , from=NORTH )
c CALL CHECKSUM( grid, APRESS, 383, "ICEDYN_DRV.f")
c CALL HALO_UPDATE(grid, APRESS, from=NORTH )
CALL CHECKSUM( grid, OGEOZA, 385, "ICEDYN_DRV.f")
CALL HALO_UPDATE(grid, OGEOZA, from=NORTH )
CALL CHECKSUM( grid, DYV , 387, "ICEDYN_DRV.f")
CALL HALO_UPDATE(grid, DYV , from=NORTH )
CALL CHECKSUM( grid, MSI , 389, "ICEDYN_DRV.f")
CALL HALO_UPDATE(grid, MSI , from=NORTH )
CALL CHECKSUM( grid, SNOWI , 391, "ICEDYN_DRV.f")
CALL HALO_UPDATE(grid, SNOWI , from=NORTH )
DO J=J_0,J_1S
DO I=1,IM
IF(FOCEAN(I,J+1).gt.0 .and. FOCEAN(I,J).gt.0. .and.
* RSI(I,J)+RSI(I,J+1).gt.0.) THEN
c PGFV(I,J)=-((APRESS(I,J+1)-APRESS(I,J))*BYRHOI
c * +OGEOZA(I,J+1)-OGEOZA(I,J))/DYV(J+1)
PGFV(I,J)=-(OGEOZA(I,J+1)-OGEOZA(I,J)+
* (RSI(I,J+1)*(MSI(I,J+1)+SNOWI(I,J+1)+ACE1I) -
* RSI(I,J)*(MSI(I,J)+SNOWI(I,J)+ACE1I))*GRAV/RHOWS )
* /DYV(J+1)
ELSE
PGFV(I,J)=0.
END IF
END DO
END DO
C**** DMUA is defined over the whole box (not just over ptype)
C**** Convert to stress over ice fraction only (on atmospheric grid)
DO J=J_0,J_1
DO I=1,IM
IF (FOCEAN(I,J)*RSI(I,J).gt.0) THEN
DMUA(I,J,2) = DMUA(I,J,2)/(FOCEAN(I,J)*RSI(I,J))
DMVA(I,J,2) = DMVA(I,J,2)/(FOCEAN(I,J)*RSI(I,J))
ELSE
DMUA(I,J,2) = 0.
DMVA(I,J,2) = 0.
END IF
END DO
END DO
C**** Set up ice grid variables
C**** HEFF,AREA on primary (tracer) grid for ice
C**** Fill halos for RSI, FOCEAN, MSI
CALL CHECKSUM( grid, RSI , 426, "ICEDYN_DRV.f")
CALL HALO_UPDATE(grid, RSI , from=NORTH+SOUTH )
CALL CHECKSUM( grid, FOCEAN, 428, "ICEDYN_DRV.f")
CALL HALO_UPDATE(grid, FOCEAN, from=NORTH+SOUTH )
CALL CHECKSUM( grid, MSI , 430, "ICEDYN_DRV.f")
CALL HALO_UPDATE(grid, MSI , from=NORTH+SOUTH )
do j=J_0NXY,J_1NXY
do i=2,NX1-1
HEFF(I,J)=RSI(I-1,J)*(ACE1I+MSI(I-1,J))*BYRHOI
AREA(I,J)=RSI(I-1,J)
c if heffm(i,j).eq.1 .and. focean(i-1,j).eq.0 interpolate.....
c if (heffm(i,j).eq.1 .and. focean(i-1,j).eq.0) then
c rsifix=(rsi(i-2,j)*focean(i-2,j)+rsi(i,j)*focean(i,j)+rsi(i-1,j+1)
c * *focean(i-1,j+1)+rsi(i-1,j-1)*focean(i-1,j-1))/(focean(i-2,j)
c * +focean(i,j)+focean(i-1,j+1)+focean(i-1,j-1))
c msifix=(msi(i-2,j)*focean(i-2,j)+msi(i,j)*focean(i,j)+msi(i-1,j+1)
c * *focean(i-1,j+1)+msi(i-1,j-1)*focean(i-1,j-1))/(focean(i-2,j)
c * +focean(i,j)+focean(i-1,j+1)+focean(i-1,j-1))
c HEFF(I,J)=rsifix*(ACE1I+msifix)*BYRHOI
c AREA(I,J)=rsifix
c end if
HEFF(I,J)=HEFF(I,J)*HEFFM(I,J)
enddo
enddo
C**** fill in overlap regions
DO J=J_0NXY,J_1NXY
HEFF(1,J)=HEFF(NX1-1,J)
AREA(1,J)=AREA(NX1-1,J)
HEFF(NX1,J)=HEFF(2,J)
AREA(NX1,J)=AREA(2,J)
END DO
C****
C**** Set up mass per unit area and coriolis term (on ice grid B)
C****
C**** Update halo for HEFF
CALL CHECKSUM(grid, HEFF, 463, "ICEDYN_DRV.f")
CALL HALO_UPDATE(grid, HEFF, from=NORTH )
DO J=J_0NXY,J_1NXYS
DO I=1,NX1-1
AMASS(I,J)=RHOI*0.25*(HEFF(I,J)
* +HEFF(I+1,J)+HEFF(I,J+1)+HEFF(I+1,J+1))
COR(I,J)=AMASS(I,J)*2.0*OMEGA*SINEN(I,J)
END DO
END DO
c**** set north pole
if (grid%HAVE_NORTH_POLE) then
do i=1,nx1
AMASS(i,jm)= RHOI*HEFF(1,JM)
COR (i,jm)= AMASS(i,jm)*2.0*OMEGA*SINEN(1,JM)
end do
end if !end NORTH_POLE block if
c**** interpolate air, current and ice velocity from C grid to B grid
C**** This should be more generally from ocean grid to ice grid
C**** NOTE: UOSURF, VOSURF are expected to be on the C-grid
C**** Update halo for USI,UOSURF,PGFU
CALL CHECKSUM(grid, USI , 486, "ICEDYN_DRV.f")
CALL HALO_UPDATE(grid, USI , from=NORTH )
CALL CHECKSUM(grid, UOSURF, 488, "ICEDYN_DRV.f")
CALL HALO_UPDATE(grid, UOSURF, from=NORTH )
CALL CHECKSUM(grid, PGFU , 490, "ICEDYN_DRV.f")
CALL HALO_UPDATE(grid, PGFU , from=NORTH )
do j=j_0,j_1s
im1=im
do i=1,im
UIB (i,j)=0.5*(USI (im1,j) +USI (im1,j+1)) ! iceC--> iceB
GWATX(i,j)=0.5*(UOSURF(im1,j)+UOSURF(im1,j+1)) ! ocnC--> iceB
PGFUB(i,j)=0.5*(PGFU(im1,j) +PGFU(im1,j+1)) ! atmC--> iceB
VIB (i,j)=0.5*(VSI (im1,j) +VSI (i,j))
GWATY(i,j)=0.5*(VOSURF(im1,j)+VOSURF(i,j))
PGFVB(i,j)=0.5*(PGFV(im1,j) +PGFV(i,j))
im1=i
enddo
enddo
c**** set north pole
if (grid%HAVE_NORTH_POLE) then
do i=1,im
UIB (i,jm)=USI(1,jm)
GWATX(i,jm)=UOSURF(1,jm)
PGFUB(i,jm)=PGFU(1,jm)
VIB (i,jm)=0.
GWATY(i,jm)=0.
PGFVB(i,jm)=0.
enddo
end if !end NORTH_POLE block if
DO J=J_0NXY,J_1NXY
UIB(nx1-1,J)=UIB(1,J)
VIB(nx1-1,J)=VIB(1,J)
GWATX(nx1-1,J)=GWATX(1,J)
GWATY(nx1-1,J)=GWATY(1,J)
PGFUB(nx1-1,J)=PGFUB(1,J)
PGFVB(nx1-1,J)=PGFVB(1,J)
UIB(nx1,J)=UIB(2,J)
VIB(nx1,J)=VIB(2,J)
GWATX(nx1,J)=GWATX(2,J)
GWATY(nx1,J)=GWATY(2,J)
PGFUB(nx1,J)=PGFUB(2,J)
PGFVB(nx1,J)=PGFVB(2,J)
END DO
c**** interpolate air stress from A grid in atmos, to B grid in ice
C**** change of unit from change of momentum, to flux
C**** Update halo for USI,UOSURF,PGFU
CALL CHECKSUM(grid, DMUA , 534, "ICEDYN_DRV.f")
CALL HALO_UPDATE(grid, DMUA , from=NORTH )
CALL CHECKSUM(grid, DMVA , 536, "ICEDYN_DRV.f")
CALL HALO_UPDATE(grid, DMVA , from=NORTH )
do j=j_0,j_1s
im1=im
do i=1,im
GAIRX(i,j)=0.25*(dmua(i,j,2)+dmua(im1,j,2)+dmua(im1,j+1,2)
& +dmua(i,j+1,2))*bydts
GAIRY(i,j)=0.25*(dmva(i,j,2)+dmva(im1,j,2)+dmva(im1,j+1,2)
& +dmva(i,j+1,2))*bydts
im1=i
enddo
enddo
IF (grid%HAVE_NORTH_POLE) THEN
GAIRX(1:nx1,jm)=dmua(1,jm,2)*bydts
GAIRY(1:nx1,jm)=dmva(1,jm,2)*bydts
END IF
do j=J_0NXY,J_1NXYS
GAIRX(nx1-1,j)=GAIRX(1,j)
GAIRY(nx1-1,j)=GAIRY(1,j)
GAIRX(nx1,j)=GAIRX(2,j)
GAIRY(nx1,j)=GAIRY(2,j)
enddo
c**** read in sea ice velocity
DO J=J_0NXY,J_1NXY
DO I=1,NX1
UICE(I,J,1)=UIB(I,J)
VICE(I,J,1)=VIB(I,J)
UICE(I,J,2)=0.
VICE(I,J,2)=0.
UICE(I,J,3)=0.
VICE(I,J,3)=0.
END DO
END DO
C**** do the looping over pseudo-timesteps
CALL VPICEDYN
C**** Calculate stress on ice velocity grid (B grid)
DO J=J_0NXY,J_1NXY
hemi=-1.
if (J.gt.NY1/2) hemi=1.
DO I=1,NX1-1
DMU(i,j)=DTS*dwatn(i,j)*(COSWAT*(UICE(i,j,1)-GWATX(i,j))-
* HEMI*SINWAT*(VICE(i,j,1)-GWATY(i,j)))
DMV(i,j)=DTS*dwatn(i,j)*(HEMI*SINWAT*(UICE(i,j,1)-GWATX(i,j))
* +COSWAT*(VICE(i,j,1)-GWATY(i,j)))
END DO
END DO
DO J=J_0NXY,J_1NXY
DMU(1,J)=DMU(NX1-1,J)
DMU(NX1,J)=DMU(2,J)
END DO
C**** interpolate ice velocity and stress from B grid to C grid in atm
C**** Update halos for UICE and DMU
CALL CHECKSUM( grid, UICE, 593, "ICEDYN_DRV.f" )
CALL HALO_UPDATE(grid, UICE, from=SOUTH )
CALL CHECKSUM( grid, DMU, 595, "ICEDYN_DRV.f" )
CALL HALO_UPDATE(grid, DMU, from=SOUTH )
do j=J_0STG,J_1STG
i=im
do ip1=1,im
usi(i,j)=0.5*(uice(i+1,j-1,1)+uice(i+1,j,1))
IF (abs(USI(I,J)).lt.1d-10) USI(I,J)=0
DMUI(I,J) = 0.5*(dmu(i+1,j-1)+dmu(i+1,j))
C**** Rescale DMUI to be net momentum into ocean
DMUI(I,J) = 0.5*DMUI(I,J)*(FOCEAN(I,J)*RSI(I,J)+FOCEAN(IP1,J)
* *RSI(IP1,J))
i=ip1
enddo
enddo
C**** Update halos for FOCEAN, DXYS, DXYV
CALL CHECKSUM( grid, FOCEAN, 610, "ICEDYN_DRV.f" )
CALL HALO_UPDATE(grid, FOCEAN, from=NORTH )
CALL CHECKSUM( grid, DXYS , 612, "ICEDYN_DRV.f" )
CALL HALO_UPDATE(grid, DXYS , from=NORTH )
CALL CHECKSUM( grid, DXYV , 614, "ICEDYN_DRV.f" )
CALL HALO_UPDATE(grid, DXYV , from=NORTH )
do j=j_0,j_1s
do i=1,im
vsi(i,j)=0.5*(vice(i,j,1)+vice(i+1,j,1))
IF (abs(VSI(I,J)).lt.1d-10) VSI(I,J)=0
DMVI(I,J) = 0.5*(dmv(i,j)+dmv(i+1,j))
C**** Rescale DMVI to be net momentum into ocean
IF (J.lt.JM-1) DMVI(I,J) = 0.5*DMVI(I,J)*(FOCEAN(I,J)*RSI(I,J)
* *DXYN(J)+FOCEAN(I,J+1)*RSI(I,J+1)*DXYS(J+1))/DXYV(J+1)
IF (J.eq.JM-1) DMVI(I,JM-1) = 0.5*DMVI(I,JM-1)*(FOCEAN(I,JM-1)
* *RSI(I,JM-1)*DXYN(JM-1)+FOCEAN(1,JM)*RSI(1,JM)*DXYS(JM))
* /DXYV(JM)
enddo
enddo
IF (grid%HAVE_SOUTH_POLE) then
usi(1:im,1)=0.
dmui(1:im,1)=0.
END IF
IF (grid%HAVE_NORTH_POLE) then
vsi(1:im,jm)=0.
dmvi(1:im,jm)=0.
END IF
C**** Calculate ustar*2*rho for ice-ocean fluxes on atmosphere grid
C**** UI2rho = | tau |
C**** Update halos for DMU, and DMV
CALL CHECKSUM( grid, DMU , 642, "ICEDYN_DRV.f" )
CALL HALO_UPDATE(grid, DMU , from=SOUTH )
CALL CHECKSUM( grid, DMV , 644, "ICEDYN_DRV.f" )
CALL HALO_UPDATE(grid, DMV , from=SOUTH )
do j=J_0,J_1
do i=1,imaxj(j)
UI2rho(i,j)=0
if (FOCEAN(I,J)*RSI(i,j).gt.0) THEN
C**** calculate 4 point average of B grid values of stresses
duA = 0.5*(DXYN(J)*(dmu(i+1,j)+dmu(i,j))+DXYS(j)*(dmu(i+1
* ,j-1)+dmu(i,j-1)))*BYDXYP(J)
dvA = 0.5*(DXYN(J)*(dmv(i+1,j)+dmv(i,j))+DXYS(j)*(dmv(i+1
* ,j-1)+dmv(i,j-1)))*BYDXYP(J)
UI2rho(i,j)= sqrt (duA**2 + dvA**2) * bydts
end if
end do
end do
C**** set north pole in C grid (atm)
IF (grid%HAVE_NORTH_POLE) THEN
USINP=0.
DMUINP=0.
do i=1,im
USINP = USINP + USI(i,jm)
DMUINP = DMUINP + DMUI(i,jm)
enddo
USINP=USINP/IM
DMUINP=DMUINP/IM
USI(1:im,jm)=USINP
DMUI(1:im,jm)=DMUINP
VSI(1:im,jm)=0.
DMVI(1:im,jm)=0.
END IF
C**** calculate mass fluxes for the ice advection
C**** Update halos for FOCEAN, and RSI
CALL CHECKSUM( grid, FOCEAN, 679, "ICEDYN_DRV.f" )
CALL HALO_UPDATE(grid, FOCEAN, from=NORTH )
CALL CHECKSUM( grid, RSI , 681, "ICEDYN_DRV.f" )
CALL HALO_UPDATE(grid, RSI , from=NORTH )
DO J=J_0,J_1S
I=IM
DO IP1=1,IM
USIDT(I,J)=0.
IF (FOCEAN(I,J).gt.0 .and. FOCEAN(IP1,J).gt.0. .and.
* RSI(I,J)+RSI(IP1,J).gt.1d-4) THEN
USIDT(I,J)=USI(I,J)*DTS
ICIJ(I,J,IJ_USI) =ICIJ(I,J,IJ_USI) +(RSI(I,J)+RSI(IP1,J))
* *USI(i,j)
ICIJ(I,J,IJ_DMUI)=ICIJ(I,J,IJ_DMUI)+DMUI(i,j)
END IF
VSIDT(I,J)=0.
IF (FOCEAN(I,J+1).gt.0 .and. FOCEAN(I,J).gt.0. .and.
* RSI(I,J)+RSI(I,J+1).gt.1d-4) THEN
VSIDT(I,J)=VSI(I,J)*DTS
ICIJ(I,J,IJ_VSI) =ICIJ(I,J,IJ_VSI) +(RSI(I,J)+RSI(I,J+1))
* *VSI(i,j)
ICIJ(I,J,IJ_DMVI)=ICIJ(I,J,IJ_DMVI)+DMVI(i,j)
END IF
ICIJ(I,J,IJ_PICE)=ICIJ(I,J,IJ_PICE)+ RSI(I,J)*press(i+1,j)
ICIJ(I,J,IJ_RSI) =ICIJ(I,J,IJ_RSI) + RSI(I,J)
I=IP1
END DO
END DO
IF (grid%HAVE_NORTH_POLE) THEN
VSIDT(1:IM,JM)=0.
USIDT(1:IM,JM)=USI(1,JM)*DTS
ICIJ(1,JM,IJ_USI) =ICIJ(1,JM,IJ_USI) +RSI(1,JM)*USI(1,JM)
ICIJ(1,JM,IJ_DMUI)=ICIJ(1,JM,IJ_DMUI)+DMUI(1,JM)
ICIJ(1,JM,IJ_RSI) =ICIJ(1,JM,IJ_RSI) +RSI(1,JM)
ICIJ(1,JM,IJ_PICE)=ICIJ(1,JM,IJ_PICE)+RSI(1,JM)*press(1,JM)
END IF
C**** Set uisurf,visurf for use in atmospheric drag calculations
uisurf=usi
visurf=vsi
C****
END SUBROUTINE DYNSI
SUBROUTINE ADVSI 1,11
!@sum ADVSI advects sea ice
!@+ Currently set up to advect ice on AGCM grid (i.e. usidt/vsidt are
!@+ on the AGCM grid, and RSI/MSI/HSI etc. are unchanged)
!@+ At some point this will change (USIDT/VSIDT on ice grid, and RSI
!@+ etc. will need to be interpolated back and forth).
!@auth Gary Russell/Gavin Schmidt
USE CONSTANT, only : byshi,lhm,grav
USE MODEL_COM, only : im,jm,focean,p,ptop,kocean
USE DOMAIN_DECOMP, only : grid, GET
USE DOMAIN_DECOMP, only : HALO_UPDATE, CHECKSUM, SOUTH, NORTH
USE GEOM, only : dxyp,dyp,dxp,dxv,bydxyp,imaxj
c USE ICEGEOM, only : dxyp,dyp,dxp,dxv,bydxyp ?????
USE ICEDYN_COM, only : usidt,vsidt,rsix,rsiy,rsisave,icij,ij_musi
* ,ij_mvsi,ij_husi,ij_hvsi,ij_susi,ij_svsi
USE SEAICE, only : ace1i,xsi
USE SEAICE_COM, only : rsi,msi,snowi,hsi,ssi,lmi
USE FLUXES, only : gtemp,apress,msicnv,fwsim
USE DAGCOM, only : oa
IMPLICIT NONE
REAL*8, DIMENSION(IM) :: FAW,FASI,FXSI,FYSI
!@var NTRICE max. number of tracers to be advected (mass/heat/salt+)
INTEGER, PARAMETER :: NTRICE=2+2*LMI
REAL*8 FMSI(NTRICE,IM),SFMSI(NTRICE),AMSI(NTRICE)
REAL*8 BYFOA(IM,grid%J_STRT_HALO:grid%J_STOP_HALO)
INTEGER I,J,L,IM1,IP1,K
REAL*8 SFASI,DMHSI,ASI,YRSI,XRSI,FRSI,SICE
!@var MHS mass/heat/salt content of sea ice
REAL*8 MHS(NTRICE,IM,grid%J_STRT_HALO:grid%J_STOP_HALO)
C****
C**** FLUXCB USIDT U compon of time integrated sea ice velocity (m)
C**** VSIDT V compon of time integrated sea ice velocity (m)
C****
C**** FAW flux of surface water area (m^2) = USIDT*DYP
C**** FASI flux of sea ice area (m^2) = USIDT*DYP*RSIedge
C**** FMSI flux of sea ice mass (kg) or heat (J) or salt (kg)
INTEGER J_0, J_1, J_0S, J_1S
C**** Get grid parameters
CALL GET(grid, J_STRT =J_0, J_STOP =J_1,
& J_STRT_SKP =J_0S, J_STOP_SKP =J_1S )
C**** Regularise ice concentration gradients to prevent advection errors
DO J=J_0S, J_1S
DO I=1,IM
IF (RSI(I,J).gt.1d-4) THEN
IF (RSISAVE(I,J).gt.RSI(I,J)) THEN ! reduce gradients
FRSI=(RSISAVE(I,J)-RSI(I,J))/RSISAVE(I,J)
RSIX(I,J)=RSIX(I,J)*(1.-FRSI)
RSIY(I,J)=RSIY(I,J)*(1.-FRSI)
END IF
IF(RSI(I,J)-RSIX(I,J).lt.0.) RSIX(I,J) = RSI(I,J)
IF(RSI(I,J)+RSIX(I,J).lt.0.) RSIX(I,J) = -RSI(I,J)
IF(RSI(I,J)-RSIX(I,J).gt.1d0) RSIX(I,J) = RSI(I,J)-1d0
IF(RSI(I,J)+RSIX(I,J).gt.1d0) RSIX(I,J) =1d0-RSI(I,J)
IF(RSI(I,J)-RSIY(I,J).lt.0.) RSIY(I,J) = RSI(I,J)
IF(RSI(I,J)+RSIY(I,J).lt.0.) RSIY(I,J) = -RSI(I,J)
IF(RSI(I,J)-RSIY(I,J).gt.1d0) RSIY(I,J) = RSI(I,J)-1d0
IF(RSI(I,J)+RSIY(I,J).gt.1d0) RSIY(I,J) =1d0-RSI(I,J)
ELSE
RSIX(I,J) = 0. ; RSIY(I,J) = 0.
END IF
END DO
END DO
C**** set up local MHS array to contain all advected quantities
C**** MHS(1:2) = MASS, MHS(3:6) = HEAT, MHS(7:10)=SALT
C**** Currently this is on atmospheric grid
MHS(1,:,:) = ACE1I + SNOWI
MHS(2,:,:) = MSI
DO L=1,LMI
MHS(L+2,:,:) = HSI(L,:,:)
MHS(L+2+LMI,:,:) = SSI(L,:,:)
END DO
C**** define inverse area array
DO J=J_0, J_1
DO I=1,IM
IF (FOCEAN(I,J).gt.0) THEN
BYFOA(I,J)=BYDXYP(J)/FOCEAN(I,J)
ELSE
BYFOA(I,J)=0.
END IF
END DO
END DO
C****
C**** North-South Advection of Sea Ice
C****
SFASI = 0.
SFMSI(1:NTRICE) = 0.
DO 340 I=1,IM
C****
C**** Calculate south-north sea ice fluxes at grid box edges
C****
C**** Update halo of DXV,RSIY,RSI,RSIX,FOCEAN,BYDXYP,and MHS
CALL CHECKSUM(grid, DXV , 854, "ICEDYN_DRV.f")
CALL HALO_UPDATE(grid, DXV , FROM=NORTH)
CALL CHECKSUM(grid, RSIY , 856, "ICEDYN_DRV.f")
CALL HALO_UPDATE(grid, RSIY , FROM=NORTH)
CALL CHECKSUM(grid, RSIX , 858, "ICEDYN_DRV.f")
CALL HALO_UPDATE(grid, RSIX , FROM=NORTH)
CALL CHECKSUM(grid, RSI , 860, "ICEDYN_DRV.f")
CALL HALO_UPDATE(grid, RSI , FROM=NORTH)
CALL CHECKSUM(grid, FOCEAN, 862, "ICEDYN_DRV.f")
CALL HALO_UPDATE(grid,FOCEAN, FROM=NORTH)
CALL CHECKSUM(grid, BYDXYP, 864, "ICEDYN_DRV.f")
CALL HALO_UPDATE(grid,BYDXYP, FROM=NORTH)
CALL CHECKSUM(grid, MHS , 866, "ICEDYN_DRV.f")
CALL HALO_UPDATE(grid, MHS , FROM=NORTH)
DO 120 J=2,JM-2
IF(VSIDT(I,J).eq.0.) GO TO 120
FAW(J) = VSIDT(I,J)*DXV(J+1) ! be careful with atm. grid index
IF(VSIDT(I,J).le.0.) THEN
C**** Sea ice velocity is southward at grid box edge
FASI(J)=FAW(J)*(RSI(I,J+1)-(1d0+FAW(J)*BYDXYP(J+1))*RSIY(I,J+1))
* *FOCEAN(I,J+1)
FXSI(J)=FAW(J)*RSIX(I,J+1)*FOCEAN(I,J+1)
FYSI(J)=FAW(J)*(FAW(J)*BYDXYP(J+1)*
* FAW(J)*RSIY(I,J+1)*FOCEAN(I,J+1) - 3d0*FASI(J))
FMSI(1:NTRICE,J) = FASI(J)*MHS(1:NTRICE,I,J+1)
ELSE
C**** Sea ice velocity is northward at grid box edge
FASI(J)=FAW(J)*(RSI(I,J)+(1d0-FAW(J)*BYDXYP(J))*RSIY(I,J))
* *FOCEAN(I,J)
FXSI(J)=FAW(J)*RSIX(I,J)*FOCEAN(I,J)
FYSI(J)=FAW(J)*(FAW(J)*BYDXYP(J)*FAW(J)*RSIY(I,J)*FOCEAN(I,J)
* -3d0*FASI(J))
FMSI(1:NTRICE,J) = FASI(J)*MHS(1:NTRICE,I,J)
END IF
ICIJ(I,J,IJ_MVSI)=ICIJ(I,J,IJ_MVSI)+SUM(FMSI(1:2,J))
ICIJ(I,J,IJ_HVSI)=ICIJ(I,J,IJ_HVSI)+SUM(FMSI(3:2+LMI,J))
ICIJ(I,J,IJ_SVSI)=ICIJ(I,J,IJ_SVSI)+SUM(FMSI(3+LMI:2+2*LMI,J))
120 CONTINUE
C****
C**** Calculate south-north sea ice fluxes near North Pole
C****
IF (grid%HAVE_NORTH_POLE) THEN
IF(VSIDT(I,JM-1).eq.0.) GO TO 200
FAW(JM-1) = VSIDT(I,JM-1)*DXV(JM) ! careful with atm.grid index!
IF(VSIDT(I,JM-1).le.0.) THEN
C**** Sea ice velocity is southward from North Pole box
FASI(JM-1) = FAW(JM-1)*RSI(1,JM)*FOCEAN(1,JM)
FXSI(JM-1) = 0.
FYSI(JM-1) = -FAW(JM-1)*FASI(JM-1)
FMSI(1:NTRICE,JM-1) = FASI(JM-1)*MHS(1:NTRICE,1,JM)
ELSE
C**** Sea ice velocity is northward into North Pole box
FASI(JM-1) = FAW(JM-1)*FOCEAN(I,JM-1)*
* (RSI(I,JM-1)+(1d0-FAW(JM-1)*BYDXYP(JM-1))*RSIY(I,JM-1))
FXSI(JM-1) = FAW(JM-1)*RSIX(I,JM-1)*FOCEAN(I,JM-1)
FYSI(JM-1) = FAW(JM-1)*(FAW(JM-1)*BYDXYP(JM-1)*
* FAW(JM-1)*RSIY(I,JM-1)*FOCEAN(I,JM-1)-3d0*FASI(JM-1))
FMSI(1:NTRICE,JM-1) = FASI(JM-1)*MHS(1:NTRICE,I,JM-1)
END IF
C**** Accumulate sea ice leaving and entering North Pole box
SFASI = SFASI + FASI(JM-1)
SFMSI(1:NTRICE) = SFMSI(1:NTRICE) + FMSI(1:NTRICE,JM-1)
ICIJ(I,JM-1,IJ_MVSI)=ICIJ(I,JM-1,IJ_MVSI)+SUM(FMSI(1:2,JM-1))
ICIJ(I,JM-1,IJ_HVSI)=ICIJ(I,JM-1,IJ_HVSI)
& +SUM(FMSI(3:2+LMI,JM-1))
ICIJ(I,JM-1,IJ_SVSI)=ICIJ(I,JM-1,IJ_SVSI)+
* SUM(FMSI(3+LMI:2+2*LMI,JM-1))
ENDIF
C****
C**** Update sea ice variables due to south-north fluxes
C****
C****Update halo of VSIDT, FASI, FAW, FOCEAN, FMSI, and FXSI
CALL CHECKSUM(grid, VSIDT, 940, "ICEDYN_DRV.f")
CALL HALO_UPDATE(grid, VSIDT, FROM=SOUTH)
CALL CHECKSUM(grid, FASI, 942, "ICEDYN_DRV.f")
CALL HALO_UPDATE(grid, FASI, FROM=SOUTH)
CALL CHECKSUM(grid, FAW, 944, "ICEDYN_DRV.f")
CALL HALO_UPDATE(grid, FAW, FROM=SOUTH)
CALL CHECKSUM(grid, FOCEAN, 946, "ICEDYN_DRV.f")
CALL HALO_UPDATE(grid, FOCEAN, FROM=SOUTH)
CALL CHECKSUM(grid, FMSI, 948, "ICEDYN_DRV.f")
CALL HALO_UPDATE(grid, FMSI, FROM=SOUTH)
CALL CHECKSUM(grid, FXSI, 950, "ICEDYN_DRV.f")
CALL HALO_UPDATE(grid, FXSI, FROM=SOUTH)
CALL CHECKSUM(grid, FYSI, 952, "ICEDYN_DRV.f")
CALL HALO_UPDATE(grid, FYSI, FROM=SOUTH)
200 DO 330 J=J_0S, J_1S
IF(VSIDT(I,J-1)) 240,210,280
C**** VSIDT(J-1)=0.
210 IF(VSIDT(I,J)) 220,330,230
C**** VSIDT(J-1)=0, VSIDT(J)<0.
220 ASI = RSI(I,J)*DXYP(J)*FOCEAN(I,J) - FASI(J)
DO 225 K=1,NTRICE
225 AMSI(K) = RSI(I,J)*DXYP(J)*MHS(K,I,J)*FOCEAN(I,J) - FMSI(K,J)
IF(ASI.gt.DXYP(J)*FOCEAN(I,J)) GO TO 320
YRSI = (RSIY(I,J)*DXYP(J)*DXYP(J)*FOCEAN(I,J) - FYSI(J)
* + 3d0*(FAW(J)*ASI-DXYP(J)*FASI(J))) / (DXYP(J)-FAW(J))
RSI(I,J) = ASI*BYFOA(I,J)
RSIY(I,J) = YRSI*BYFOA(I,J)
RSIX(I,J) = RSIX(I,J) - FXSI(J)*BYFOA(I,J)
IF (ASI.gt.0) MHS(1:NTRICE,I,J) = AMSI(1:NTRICE)/ASI
GO TO 310
C**** VSIDT(J-1)=0, VSIDT(J)>0.
230 RSI(I,J) = RSI(I,J) - FASI(J)*BYFOA(I,J)
RSIX(I,J) = RSIX(I,J)*(1d0-FAW(J)*BYDXYP(J))
RSIY(I,J) = RSIY(I,J)*(1d0-FAW(J)*BYDXYP(J))**2
GO TO 310
C**** VSIDT(J-1)<0.
240 IF(VSIDT(I,J)) 260,250,270
C**** VSIDT(J-1)<0, VSIDT(J)=0.
250 RSI(I,J) = RSI(I,J) + FASI(J-1)*BYFOA(I,J)
RSIX(I,J) = RSIX(I,J)*(1d0+FAW(J-1)*FOCEAN(I,J-1)*BYFOA(I,J))
RSIY(I,J) = RSIY(I,J)*(1d0+FAW(J-1)*FOCEAN(I,J-1)*BYFOA(I,J))**2
GO TO 310
C**** VSIDT(J-1)<0, VSIDT(J)<0 or VSIDT(J-1)>0, VSIDT(J) not 0.
260 ASI = RSI(I,J)*DXYP(J)*FOCEAN(I,J) + ( FASI(J-1)- FASI(J))
DO 265 K=1,NTRICE
265 AMSI(K) = RSI(I,J)*DXYP(J)*MHS(K,I,J)*FOCEAN(I,J) +
* (FMSI(K,J-1)-FMSI(K,J))
IF(ASI.gt.DXYP(J)*FOCEAN(I,J)) GO TO 320
YRSI = (RSIY(I,J)*DXYP(J)*DXYP(J)*FOCEAN(I,J)+(FYSI(J-1)-FYSI(J))
* + 3d0*((FAW(J-1)+FAW(J))*ASI-DXYP(J)*(FASI(J-1)+FASI(J))))
* / (DXYP(J) + (FAW(J-1)-FAW(J)))
RSI(I,J) = ASI*BYFOA(I,J)
RSIY(I,J) = YRSI*BYFOA(I,J)
RSIX(I,J) = RSIX(I,J) + (FXSI(J-1)-FXSI(J))*BYFOA(I,J)
IF (ASI.gt.0) MHS(1:NTRICE,I,J) = AMSI(1:NTRICE)/ASI
GO TO 310
C**** VSIDT(J-1)<0, VSIDT(J)>0.
270 RSI(I,J) = RSI(I,J) + (FASI(J-1)-FASI(J))*BYFOA(I,J)
RSIX(I,J) = RSIX(I,J)*(1d0+(FAW(J-1)*FOCEAN(I,J-1)
* -FAW(J)*FOCEAN(I,J))*BYFOA(I,J))
RSIY(I,J) = RSIY(I,J)*(1d0+(FAW(J-1)*FOCEAN(I,J-1)
* -FAW(J)*FOCEAN(I,J))*BYFOA(I,J))**2
GO TO 310
C**** VSIDT(J-1)>0.
280 IF(VSIDT(I,J).ne.0.) GO TO 260
C**** VSIDT(J-1)>0, VSIDT(J)=0.
ASI = RSI(I,J)*DXYP(J)*FOCEAN(I,J) + FASI(J-1)
DO 285 K=1,NTRICE
285 AMSI(K) = RSI(I,J)*DXYP(J)*MHS(K,I,J)*FOCEAN(I,J) + FMSI(K,J-1)
IF(ASI.gt.DXYP(J)*FOCEAN(I,J)) GO TO 320
YRSI = (RSIY(I,J)*DXYP(J)*DXYP(J)*FOCEAN(I,J) + FYSI(J-1)
* + 3d0*(FAW(J-1)*ASI-DXYP(J)*FASI(J-1))) / (DXYP(J)+FAW(J-1))
RSI(I,J) = ASI*BYFOA(I,J)
RSIY(I,J) = YRSI*BYFOA(I,J)
RSIX(I,J) = RSIX(I,J) + FXSI(J-1)*BYFOA(I,J)
IF (ASI.gt.0) MHS(1:NTRICE,I,J) = AMSI(1:NTRICE)/ASI
C**** Limit RSIX and RSIY so that sea ice is positive at the edges
310 RSI(I,J) = MAX(0d0,RSI(I,J))
IF(RSI(I,J)-RSIX(I,J).lt.0.) RSIX(I,J) = RSI(I,J)
IF(RSI(I,J)+RSIX(I,J).lt.0.) RSIX(I,J) = -RSI(I,J)
IF(RSI(I,J)-RSIX(I,J).gt.1d0) RSIX(I,J) = RSI(I,J)-1d0
IF(RSI(I,J)+RSIX(I,J).gt.1d0) RSIX(I,J) =1d0-RSI(I,J)
IF(RSI(I,J)-RSIY(I,J).lt.0.) RSIY(I,J) = RSI(I,J)
IF(RSI(I,J)+RSIY(I,J).lt.0.) RSIY(I,J) = -RSI(I,J)
IF(RSI(I,J)-RSIY(I,J).gt.1d0) RSIY(I,J) = RSI(I,J)-1d0
IF(RSI(I,J)+RSIY(I,J).gt.1d0) RSIY(I,J) =1d0-RSI(I,J)
GO TO 330
C**** Sea ice crunches into itself and completely covers grid box
320 RSI(I,J) = 1d0
RSIX(I,J) = 0.
RSIY(I,J) = 0.
MHS(1,I,J) = AMSI(1)/ASI
MHS(2,I,J) =(AMSI(1)+AMSI(2))*BYFOA(I,J) - MHS(1,I,J)
DO K=1,(NTRICE-2)/LMI
MHS(3+LMI*(K-1),I,J) = AMSI(3+LMI*(K-1)) / ASI
MHS(4+LMI*(K-1),I,J) = AMSI(4+LMI*(K-1)) / ASI
DMHSI = (AMSI(3+LMI*(K-1))+AMSI(4+LMI*(K-1))+AMSI(5+LMI*(K-1))
* +AMSI(6+LMI*(K-1)))*(BYFOA(I,J) -1d0 / ASI )
MHS(5+LMI*(K-1),I,J) = AMSI(5+LMI*(K-1)) / ASI +
* XSI(3)*DMHSI
MHS(6+LMI*(K-1),I,J) = AMSI(6+LMI*(K-1)) / ASI +
* XSI(4)*DMHSI
END DO
C**** End of loop over J
330 CONTINUE
C**** End of loop over I
340 CONTINUE
C****
C**** Advection of Sea Ice leaving and entering North Pole box
C****
IF (grid%HAVE_NORTH_POLE) THEN
ASI = RSI(1,JM)*DXYP(JM)*FOCEAN(1,JM) + SFASI/IM
DO 345 K=1,NTRICE
345 AMSI(K) = RSI(1,JM)*DXYP(JM)*MHS(K,1,JM)*FOCEAN(1,JM)+
& SFMSI(K)/IM
IF(ASI.gt.DXYP(JM)*FOCEAN(1,JM)) GO TO 350
RSI(1,JM) = ASI*BYFOA(1,JM)
IF (ASI.gt.0) MHS(1:NTRICE,1,JM) = AMSI(1:NTRICE)/ASI
GO TO 400
C**** Sea ice crunches into itself at North Pole box
350 RSI(1,JM) = 1d0
MHS(1,1,JM) = AMSI(1)/ASI
MHS(2,1,JM) =(AMSI(1)+AMSI(2))*BYFOA(1,JM)-MHS(1,1,JM)
DO K=1,(NTRICE-2)/LMI
MHS(3+LMI*(K-1),1,JM) = AMSI(3+LMI*(K-1)) / ASI
MHS(4+LMI*(K-1),1,JM) = AMSI(4+LMI*(K-1)) / ASI
DMHSI = (AMSI(3+LMI*(K-1))+AMSI(4+LMI*(K-1))+
& AMSI(5+LMI*(K-1))
* +AMSI(6+LMI*(K-1)))*(BYFOA(1,JM) -1d0/ ASI)
MHS(5+LMI*(K-1),1,JM) = AMSI(5+LMI*(K-1)) / ASI +
* XSI(3)*DMHSI
MHS(6+LMI*(K-1),1,JM) = AMSI(6+LMI*(K-1)) / ASI +
* XSI(4)*DMHSI
END DO
END IF !HAVE_NORTH_POLE
C****
C**** East-West Advection of Sea Ice
C****
400 DO 640 J=J_0S, J_1S
C****
C**** Calculate west-east sea ice fluxes at grid box edges
C****
I=IM
DO IP1=1,IM
IF(USIDT(I,J).eq.0.) GO TO 420
FAW(I) = USIDT(I,J)*DYP(J)
IF(USIDT(I,J).le.0.) THEN
C**** Sea ice velocity is westward at grid box edge
FASI(I)=FAW(I)*(RSI(IP1,J)-(1d0+FAW(I)*BYDXYP(J))*RSIX(IP1,J))
* *FOCEAN(IP1,J)
FXSI(I)=FAW(I)*(FAW(I)*BYDXYP(J)*
* FAW(I)*RSIX(IP1,J)*FOCEAN(IP1,J)-3d0*FASI(I))
FYSI(I)=FAW(I)*RSIY(IP1,J)*FOCEAN(IP1,J)
FMSI(1:NTRICE,I) = FASI(I)*MHS(1:NTRICE,IP1,J)
ELSE
C**** Sea ice velocity is eastward at grid box edge
FASI(I)=FAW(I)*(RSI(I,J)+(1d0-FAW(I)*BYDXYP(J))*RSIX(I,J))
* *FOCEAN(I,J)
FXSI(I)=FAW(I)*(FAW(I)*BYDXYP(J)*FAW(I)*RSIX(I,J)*FOCEAN(I,J)
* -3d0*FASI(I))
FYSI(I)=FAW(I)*RSIY(I,J)*FOCEAN(I,J)
FMSI(1:NTRICE,I) = FASI(I)*MHS(1:NTRICE,I,J)
END IF
ICIJ(I,J,IJ_MUSI)=ICIJ(I,J,IJ_MUSI)+SUM(FMSI(1:2,I))
ICIJ(I,J,IJ_HUSI)=ICIJ(I,J,IJ_HUSI)+SUM(FMSI(3:2+LMI,I))
ICIJ(I,J,IJ_SUSI)=ICIJ(I,J,IJ_SUSI)+SUM(FMSI(3+LMI:2+2*LMI,I))
420 I=IP1
END DO
C****
C**** Update sea ice variables due to west-east fluxes
C****
IM1=IM
DO 630 I=1,IM
IF(USIDT(IM1,J)) 540,510,580
C**** USIDT(IM1)=0.
510 IF(USIDT(I,J)) 520,630,530
C**** USIDT(IM1)=0, USIDT(I)<0.
520 ASI = RSI(I,J)*DXYP(J)*FOCEAN(I,J) - FASI(I)
DO 525 K=1,NTRICE
525 AMSI(K) = RSI(I,J)*DXYP(J)*MHS(K,I,J)*FOCEAN(I,J) - FMSI(K,I)
IF(ASI.gt.DXYP(J)*FOCEAN(I,J)) GO TO 620
XRSI = (RSIX(I,J)*DXYP(J)*DXYP(J)*FOCEAN(I,J) - FXSI(I)
* + 3d0*(FAW(I)*ASI-DXYP(J)*FASI(I))) / (DXYP(J)-FAW(I))
RSI(I,J) = ASI*BYFOA(I,J)
RSIX(I,J) = XRSI*BYFOA(I,J)
RSIY(I,J) = RSIY(I,J) - FYSI(I)*BYFOA(I,J)
IF (ASI.gt.0) MHS(1:NTRICE,I,J) = AMSI(1:NTRICE)/ASI
GO TO 610
C**** USIDT(IM1)=0, USIDT(I)>0.
530 RSI(I,J) = RSI(I,J) - FASI(I)*BYFOA(I,J)
RSIX(I,J) = RSIX(I,J)*(1d0-FAW(I)*BYDXYP(J))**2
RSIY(I,J) = RSIY(I,J)*(1d0-FAW(I)*BYDXYP(J))
GO TO 610
C**** USIDT(IM1)<0.
540 IF(USIDT(I,J)) 560,550,570
C**** USIDT(IM1)<0, USIDT(I)=0.
550 RSI(I,J) = RSI(I,J) + FASI(IM1)*BYFOA(I,J)
RSIX(I,J) = RSIX(I,J)*(1d0+FAW(IM1)*FOCEAN(IM1,J)*BYFOA(I,J))**2
RSIY(I,J) = RSIY(I,J)*(1d0+FAW(IM1)*FOCEAN(IM1,J)*BYFOA(I,J))
GO TO 610
C**** USIDT(IM1)<0, USIDT(I)<0 or USIDT(IM1)>0, USIDT(I) not 0.
560 ASI = RSI(I,J)*DXYP(J)*FOCEAN(I,J) + (FASI(IM1)- FASI(I))
DO 565 K=1,NTRICE
565 AMSI(K) = RSI(I,J)*DXYP(J)*MHS(K,I,J)*FOCEAN(I,J) +
* (FMSI(K,IM1)-FMSI(K,I))
IF(ASI.gt.DXYP(J)*FOCEAN(I,J)) GO TO 620
XRSI = (RSIX(I,J)*DXYP(J)*DXYP(J)*FOCEAN(I,J)+(FXSI(IM1)-FXSI(I))
* + 3d0*((FAW(IM1)+FAW(I))*ASI-DXYP(J)*(FASI(IM1)+FASI(I))))
* / (DXYP(J) + (FAW(IM1)-FAW(I)))
RSI(I,J) = ASI*BYFOA(I,J)
RSIX(I,J) = XRSI*BYFOA(I,J)
RSIY(I,J) = RSIY(I,J) + (FYSI(IM1)-FYSI(I))*BYFOA(I,J)
IF (ASI.gt.0) MHS(1:NTRICE,I,J) = AMSI(1:NTRICE)/ASI
GO TO 610
C**** USIDT(IM1)<0, USIDT(I)>0.
570 RSI(I,J) = RSI(I,J) + (FASI(IM1)-FASI(I))*BYFOA(I,J)
RSIX(I,J) = RSIX(I,J)*(1d0+(FAW(IM1)*FOCEAN(IM1,J)
* -FAW(I)*FOCEAN(I,J))*BYFOA(I,J))**2
RSIY(I,J) = RSIY(I,J)*(1d0+(FAW(IM1)*FOCEAN(IM1,J)
* -FAW(I)*FOCEAN(I,J))*BYFOA(I,J))
GO TO 610
C**** USIDT(IM1)>0.
580 IF(USIDT(I,J).ne.0.) GO TO 560
C**** USIDT(IM1)>0, USIDT(I)=0.
ASI = RSI(I,J)*DXYP(J)*FOCEAN(I,J) + FASI(IM1)
DO 585 K=1,NTRICE
585 AMSI(K) = RSI(I,J)*DXYP(J)*MHS(K,I,J)*FOCEAN(I,J) + FMSI(K,IM1)
IF(ASI.gt.DXYP(J)*FOCEAN(I,J)) GO TO 620
XRSI = (RSIX(I,J)*DXYP(J)*DXYP(J)*FOCEAN(I,J) + FXSI(IM1)
* + 3d0*(FAW(IM1)*ASI-DXYP(J)*FASI(IM1))) / (DXYP(J)+FAW(IM1))
RSI(I,J) = ASI*BYFOA(I,J)
RSIX(I,J) = XRSI*BYFOA(I,J)
RSIY(I,J) = RSIY(I,J) + FYSI(IM1)*BYFOA(I,J)
IF (ASI.gt.0) MHS(1:NTRICE,I,J) = AMSI(1:NTRICE)/ASI
C**** Limit RSIX and RSIY so that sea ice is positive at the edges
610 RSI(I,J) = MAX(0d0,RSI(I,J))
IF(RSI(I,J)-RSIX(I,J).lt.0.) RSIX(I,J) = RSI(I,J)
IF(RSI(I,J)+RSIX(I,J).lt.0.) RSIX(I,J) = -RSI(I,J)
IF(RSI(I,J)-RSIX(I,J).gt.1d0) RSIX(I,J) = RSI(I,J)-1d0
IF(RSI(I,J)+RSIX(I,J).gt.1d0) RSIX(I,J) =1d0-RSI(I,J)
IF(RSI(I,J)-RSIY(I,J).lt.0.) RSIY(I,J) = RSI(I,J)
IF(RSI(I,J)+RSIY(I,J).lt.0.) RSIY(I,J) = -RSI(I,J)
IF(RSI(I,J)-RSIY(I,J).gt.1d0) RSIY(I,J) = RSI(I,J)-1d0
IF(RSI(I,J)+RSIY(I,J).gt.1d0) RSIY(I,J) =1d0-RSI(I,J)
GO TO 630
C**** Sea ice crunches into itself and completely covers grid box
620 RSI(I,J) = 1d0
RSIX(I,J) = 0.
RSIY(I,J) = 0.
MHS(1,I,J) = AMSI(1)/ASI
MHS(2,I,J) =(AMSI(1)+AMSI(2))*BYFOA(I,J) - MHS(1,I,J)
DO K=1,(NTRICE-2)/LMI
MHS(3+LMI*(K-1),I,J) = AMSI(3+LMI*(K-1)) / ASI
MHS(4+LMI*(K-1),I,J) = AMSI(4+LMI*(K-1)) / ASI
DMHSI = (AMSI(3+LMI*(K-1))+AMSI(4+LMI*(K-1))+AMSI(5+LMI*(K-1))
* +AMSI(6+LMI*(K-1)))*(BYFOA(I,J) -1d0/ ASI)
MHS(5+LMI*(K-1),I,J) = AMSI(5+LMI*(K-1)) / ASI +
* XSI(3)*DMHSI
MHS(6+LMI*(K-1),I,J) = AMSI(6+LMI*(K-1)) / ASI +
* XSI(4)*DMHSI
END DO
C**** End of loop over I
630 IM1=I
C**** End of loop over J
640 CONTINUE
IF (KOCEAN.eq.1) THEN ! full ocean calculation, adjust sea ice
C**** set global variables from local array
C**** Currently on atmospheric grid, so no interpolation necessary
DO J=J_0, J_1
DO I=1,IMAXJ(J)
C**** Fresh water sea ice mass convergence (needed for qflux model)
MSICNV(I,J) = RSI(I,J)*(MHS(1,I,J)+MHS(2,I,J)-SUM(MHS(3
* +LMI:2*LMI+2,I,J))) - RSISAVE(I,J)*(ACE1I+SNOWI(I,J)
* +MSI(I,J)-SUM(SSI(:,I,J)))
C**** sea ice prognostic variables
SNOWI(I,J)= MAX(0d0,MHS(1,I,J) - ACE1I)
MSI(I,J) = MHS(2,I,J)
DO L=1,LMI
HSI(L,I,J) = MHS(L+2,I,J)
END DO
C**** ensure that salinity is only associated with ice
SICE=MHS(1+2+LMI,I,J)+MHS(2+2+LMI,I,J)
IF (SNOWI(I,J).gt.XSI(2)*(ACE1I+SNOWI(I,J))) THEN
SSI(1,I,J)=0.
ELSE
SSI(1,I,J)=SICE*(XSI(1)*ACE1I-XSI(2)*SNOWI(I,J))/ACE1I
END IF
SSI(2,I,J)=SICE-SSI(1,I,J)
DO L=3,LMI
SSI(L,I,J) = MHS(L+2+LMI,I,J)
END DO
FWSIM(I,J)=RSI(I,J)*(ACE1I+SNOWI(I,J)+MSI(I,J)-
* SUM(SSI(1:LMI,I,J)))
END DO
END DO
C**** Set atmospheric arrays
DO J=J_0, J_1
DO I=1,IMAXJ(J)
IF (FOCEAN(I,J).gt.0) THEN
C**** set total atmopsheric pressure anomaly in case needed by ocean
APRESS(I,J) = 100.*( * *(SNOWI(I,J)+ACE1I+MSI(I,J))*GRAV
GTEMP(1:2,2,I,J)=((HSI(1:2,I,J)-SSI(1:2,I,J)*LHM)/
* (XSI(1:2)*(SNOWI(I,J)+ACE1I))+LHM)*BYSHI
END IF
END DO
END DO
DO I=2,IM ! North pole
APRESS(I,JM)=APRESS(1,JM)
GTEMP(1:2,2,I,JM)= GTEMP(1:2,2,1,JM)
END DO
ELSE ! fixed SST case, save implied heat convergence
DO J=J_0, J_1
DO I=1,IMAXJ(J)
IF (FOCEAN(I,J).gt.0) THEN
OA(I,J,13)=OA(I,J,13)+(RSI(I,J)*SUM(MHS(3:2+LMI,I,J))
* -RSISAVE(I,J)*SUM(HSI(1:LMI,I,J)))
C**** reset sea ice concentration
RSI(I,J)=RSISAVE(I,J)
END IF
END DO
END DO
END IF
C****
RETURN
END SUBROUTINE ADVSI
SUBROUTINE init_icedyn(iniOCEAN) 1,8
!@sum init_icedyn initializes ice dynamics variables
!@auth Gavin Schmidt
USE MODEL_COM, only : im,jm,dtsrc,foceanA=>focean
USE DAGCOM, only : ia_src
USE ICEDYN_COM
USE ICEDYN, only : setup_icedyn_grid,focean,osurf_tilt
USE FLUXES, only : uisurf,visurf
USE PARAM
IMPLICIT NONE
LOGICAL, INTENT(IN) :: iniOCEAN
INTEGER k,i,j
C**** setup ice dynamics grid
C**** Currently using ATM grid:
C**** -------------------------
C**** If a different grid (but still lat/lon) is required, edit
C**** definition of focean here and the definition of imic,jmic in
C**** ICEDYN.f, and obviously the code in DYNSI.
C**** -------------------------
C**** Set land masks for ice dynamics
focean(:,:)=foceanA(:,:) ! EDIT FOR GRID CHANGE
call setup_icedyn_grid
C**** Initiallise ice dynamics if ocean model needs initialising
if (iniOCEAN) THEN
RSIX=0.
RSIY=0.
USI=0.
VSI=0.
end if
C**** set uisurf,visurf for atmopsherice drag calculations
uisurf=usi
visurf=vsi
C**** set properties for ICIJ diagnostics
k=0
k=k+1
IJ_USI=k
lname_icij(k)="Sea ice EW velocity x POICEU"
sname_icij(k)="icij_usi"
units_icij(k)="m/s"
ia_icij(k)=ia_src
scale_icij(k)=1.
ijgrid_icij(k)=2
k=k+1
IJ_VSI=k
lname_icij(k)="Sea ice NS velocity x POICEV"
sname_icij(k)="icij_vsi"
units_icij(k)="m/s"
ia_icij(k)=ia_src
scale_icij(k)=1.
ijgrid_icij(k)=2
k=k+1
IJ_DMUI=k
lname_icij(k)="Ice-ocean EW stress"
sname_icij(k)="icij_dmui"
units_icij(k)="kg/m s^2"
ia_icij(k)=ia_src
scale_icij(k)=1./dtsrc
ijgrid_icij(k)=1
k=k+1
IJ_DMVI=k
lname_icij(k)="Ice-ocean NS stress"
sname_icij(k)="icij_dmvi"
units_icij(k)="kg/m s^2"
ia_icij(k)=ia_src
scale_icij(k)=1./dtsrc
ijgrid_icij(k)=1
k=k+1
IJ_PICE=k
lname_icij(k)="Sea ice internal pressure x POICE"
sname_icij(k)="icij_psi"
units_icij(k)="10^3 kg/m s^2"
ia_icij(k)=ia_src
scale_icij(k)=1d-3
ijgrid_icij(k)=1
k=k+1
IJ_MUSI=k
lname_icij(k)="Sea ice NS mass flux"
sname_icij(k)="icij_musi"
units_icij(k)="10^7 kg/s"
ia_icij(k)=ia_src
scale_icij(k)=1d-7/dtsrc
ijgrid_icij(k)=2
k=k+1
IJ_MVSI=k
lname_icij(k)="Sea ice EW mass flux"
sname_icij(k)="icij_mvsi"
units_icij(k)="10^7 kg/s"
ia_icij(k)=ia_src
scale_icij(k)=1d-7/dtsrc
ijgrid_icij(k)=2
k=k+1
IJ_HUSI=k
lname_icij(k)="Sea ice NS heat flux"
sname_icij(k)="icij_husi"
units_icij(k)="10^12 W"
ia_icij(k)=ia_src
scale_icij(k)=1d-12/dtsrc
ijgrid_icij(k)=2
k=k+1
IJ_HVSI=k
lname_icij(k)="Sea ice EW heat flux"
sname_icij(k)="icij_hvsi"
units_icij(k)="10^12 W"
ia_icij(k)=ia_src
scale_icij(k)=1d-12/dtsrc
ijgrid_icij(k)=2
k=k+1
IJ_SUSI=k
lname_icij(k)="Sea ice NS salt flux"
sname_icij(k)="icij_susi"
units_icij(k)="10^3 kg/s"
ia_icij(k)=ia_src
scale_icij(k)=1d-3/dtsrc
ijgrid_icij(k)=2
k=k+1
IJ_SVSI=k
lname_icij(k)="Sea ice EW salt flux"
sname_icij(k)="icij_svsi"
units_icij(k)="10^3 kg/s"
ia_icij(k)=ia_src
scale_icij(k)=1d-3/dtsrc
ijgrid_icij(k)=2
k=k+1
IJ_RSI=k
lname_icij(k)="Ocean ice fraction (ice dynamic grid)"
sname_icij(k)="icij_rsi"
units_icij(k)="%"
ia_icij(k)=ia_src
scale_icij(k)=100.
ijgrid_icij(k)=1
if (k.gt.KICIJ) then
write(6,*) "Too many ICIJ diags: increase KICIJ to at least"
* ,k
call stop_model("ICIJ diagnostic error",255)
end if
RETURN
END SUBROUTINE init_icedyn
SUBROUTINE diag_ICEDYN 1,9
!@sum diag_ICEDYN prints out diagnostics for ice dynamics
!@auth Gavin Schmidt
USE CONSTANT, only : undef,teeny
USE MODEL_COM, only : xlabel,lrunid,jmon0,jyear0,idacc,jdate0
* ,amon0,jdate,amon,jyear
USE DAGCOM, only : qdiag,acc_period
USE ICEDYN_COM
USE ICEDYN, only : focean
USE FILEMANAGER, only : openunit
IMPLICIT NONE
REAL*8, DIMENSION(IMIC,JMIC) :: Q,ADENOM
INTEGER I,J,L,N,KXLB,ijgrid,IP1,k1,k
CHARACTER XLB*30
CHARACTER TITLE*80,lname*50,sname*30,units*50
character*50 :: unit_string
REAL*8 QJ(JM),QSUM
REAL*8 byiacc
IF (.not. QDIAG) RETURN
C**** determine label to be added to all titles
KXLB = INDEX(XLABEL(1:11),'(')-1
IF(KXLB.le.0) KXLB = 10
XLB = ' '
XLB(1:13)=acc_period(1:3)//' '//acc_period(4:12)
XLB = TRIM(XLB)//" "//XLABEL(1:KXLB)
C**** Open output files
call open_ij(trim(acc_period)//'.icij'//
* XLABEL(1:LRUNID),imic,jmic)
C****
C**** Simple scaled ICIJ diagnostics
C****
DO K=1,KICIJ
byiacc=1./(IDACC(IA_ICIJ(K))+teeny)
adenom=1.
lname=lname_icij(k)
k1 = index(lname,' x ')
if (k1 .gt. 0) then
if (index(lname,' x POICE ') .gt. 0) then
do j=1,jmic
do i=1,imic
adenom(i,j)=icij(i,j,ij_rsi) * byiacc
end do
end do
else if (index(lname,' x POICEU') .gt. 0) then
do j=1,jmic
i=imic
do ip1=1,imic
adenom(i,j)=0.5*(icij(i,j,ij_rsi)+icij(ip1,j,ij_rsi))
* * byiacc
i=ip1
end do
end do
else if (index(lname,' x POICEV') .gt. 0) then
do j=1,jmic-1
do i=1,imic
adenom(i,j)=0.5*(icij(i,j,ij_rsi)+icij(i,j+1,ij_rsi))
* * byiacc
end do
end do
end if
lname(k1:50) = ' '
end if
Q=UNDEF ; QJ=UNDEF ; QSUM=UNDEF
DO J=1,JMIC
DO I=1,IMIC
IF (ADENOM(I,J).gt.0 .and. FOCEAN(I,J).gt.0.5)
* Q(I,J)=SCALE_ICIJ(K)*ICIJ(I,J,K)*byiacc/adenom(i,j)
END DO
END DO
Q(2:IMIC,JMIC)=Q(1,JMIC)
Q(2:IMIC,1)=Q(1,1)
TITLE=trim(LNAME)//" ("//trim(UNITS_ICIJ(K))//") "
TITLE(51:80)=XLB
CALL POUT_IJ(TITLE,SNAME_ICIJ(K),LNAME_ICIJ(K),UNITS_ICIJ(K),Q
* ,QJ,QSUM,IJGRID_ICIJ(K),IJGRID_ICIJ(K)) ! assume igrid=jgrid
END DO
call close_ij
C****
RETURN
END SUBROUTINE diag_ICEDYN