C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/seaice/seaice_init_varia.F,v 1.7 2007/06/05 16:01:22 mlosch Exp $ C $Name: $ #include "SEAICE_OPTIONS.h" CStartOfInterface SUBROUTINE SEAICE_INIT_VARIA( myThid ) C /==========================================================\ C | SUBROUTINE SEAICE_INIT_VARIA | C | o Initialization of sea ice model. | C |==========================================================| C \==========================================================/ IMPLICIT NONE C === Global variables === #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "GRID.h" #include "SEAICE.h" CML#include "SEAICE_GRID.h" #include "SEAICE_DIAGS.h" #include "SEAICE_PARAMS.h" #include "FFIELDS.h" #ifdef ALLOW_EXCH2 #include "W2_EXCH2_TOPOLOGY.h" #include "W2_EXCH2_PARAMS.h" #endif /* ALLOW_EXCH2 */ C === Routine arguments === C myThid - Thread no. that called this routine. INTEGER myThid CEndOfInterface C === Local variables === C i,j,k,bi,bj - Loop counters INTEGER i, j, k, bi, bj _RS mask_uice INTEGER myIter, myTile cph( cph make sure TAF sees proper initialisation cph to avoid partial recomputation issues DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) c DO K=1,3 DO J=1-OLy,sNy+OLy DO I=1-OLx,sNx+OLx HEFF(I,J,k,bi,bj)=ZERO AREA(I,J,k,bi,bj)=ZERO UICE(I,J,k,bi,bj)=ZERO VICE(I,J,k,bi,bj)=ZERO ENDDO ENDDO ENDDO c DO J=1-OLy,sNy+OLy DO I=1-OLx,sNx+OLx HSNOW(I,J,bi,bj)=ZERO ZETA(I,J,bi,bj)=ZERO ENDDO ENDDO c ENDDO ENDDO cph) C-- Initialize grid info DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx HEFFM(i,j,bi,bj)=ONE IF (_hFacC(i,j,1,bi,bj).eq.0.) HEFFM(i,j,bi,bj)=ZERO ENDDO ENDDO DO J=1-OLy+1,sNy+OLy DO I=1-OLx+1,sNx+OLx #ifndef SEAICE_CGRID UVM(i,j,bi,bj)=ZERO mask_uice=HEFFM(I,J, bi,bj)+HEFFM(I-1,J-1,bi,bj) & +HEFFM(I,J-1,bi,bj)+HEFFM(I-1,J, bi,bj) IF(mask_uice.GT.3.5) UVM(I,J,bi,bj)=ONE #else seaiceMaskU(I,J,bi,bj)= 0.0 _d 0 seaiceMaskV(I,J,bi,bj)= 0.0 _d 0 mask_uice=HEFFM(I,J,bi,bj)+HEFFM(I-1,J ,bi,bj) IF(mask_uice.GT.1.5) seaiceMaskU(I,J,bi,bj)=ONE mask_uice=HEFFM(I,J,bi,bj)+HEFFM(I ,J-1,bi,bj) IF(mask_uice.GT.1.5) seaiceMaskV(I,J,bi,bj)=ONE #endif /* not SEAICE_CGRID */ ENDDO ENDDO #ifdef ALLOW_EXCH2 #ifdef SEAICE_CGRID #else C-- Special stuff for cubed sphere: assume grid is rectangular and C set UV mask to zero except for Arctic and Antarctic cube faces. IF (useCubedSphereExchange) THEN myTile = W2_myTileList(bi) IF ( exch2_myFace(myTile) .EQ. 1 .OR. & exch2_myFace(myTile) .EQ. 2 .OR. & exch2_myFace(myTile) .EQ. 4 .OR. & exch2_myFace(myTile) .EQ. 5 ) THEN DO J=1-OLy,sNy+OLy DO I=1-OLx,sNx+OLx UVM(i,j,bi,bj)=ZERO ENDDO ENDDO ELSEIF ( exch2_isWedge(myTile) .EQ. 1 ) THEN I=1 DO J=1-OLy,sNy+OLy UVM(i,j,bi,bj)=ZERO ENDDO ELSEIF ( exch2_isSedge(myTile) .EQ. 1 ) THEN J=1 DO I=1-OLx,sNx+OLx UVM(i,j,bi,bj)=ZERO ENDDO ENDIF ENDIF #endif #endif /* ALLOW_EXCH2 */ DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx TICE(I,J,bi,bj)=273.0 _d 0 #ifdef SEAICE_MULTICATEGORY DO k=1,MULTDIM TICES(I,J,k,bi,bj)=273.0 _d 0 ENDDO #endif /* SEAICE_MULTICATEGORY */ UICEC (I,J,bi,bj)=ZERO VICEC (I,J,bi,bj)=ZERO DWATN (I,J,bi,bj)=ZERO #ifndef SEAICE_CGRID DAIRN (I,J,bi,bj)=ZERO AMASS (I,J,bi,bj)=1000.0 _d 0 #else seaiceMassC(I,J,bi,bj)=1000.0 _d 0 seaiceMassU(I,J,bi,bj)=1000.0 _d 0 seaiceMassV(I,J,bi,bj)=1000.0 _d 0 #endif GWATX (I,J,bi,bj)=ZERO GWATY (I,J,bi,bj)=ZERO ENDDO ENDDO C-- Choose a proxy level for geostrophic velocity, DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx #ifdef SEAICE_TEST_ICE_STRESS_1 KGEO(I,J,bi,bj) = 1 #else /* SEAICE_TEST_ICE_STRESS_1 */ IF (klowc(i,j,bi,bj) .LT. 2) THEN KGEO(I,J,bi,bj) = 1 ELSE KGEO(I,J,bi,bj) = 2 DO WHILE ( abs(rC(KGEO(I,J,bi,bj))) .LT. 50.0 .AND. & KGEO(I,J,bi,bj) .LT. (klowc(i,j,bi,bj)-1) ) KGEO(I,J,bi,bj) = KGEO(I,J,bi,bj) + 1 ENDDO ENDIF #endif /* SEAICE_TEST_ICE_STRESS_1 */ ENDDO ENDDO ENDDO ENDDO C-- Update overlap regions #ifdef SEAICE_CGRID CALL EXCH_UV_XY_RL(seaiceMaskU,seaiceMaskV,.FALSE.,myThid) #else _EXCH_XY_R8(UVM, myThid) #endif C-- Now lets look at all these beasts IF ( debugLevel .GE. debLevB ) THEN myIter=0 CALL PLOT_FIELD_XYRL( HEFFM , 'Current HEFFM ' , & myIter, myThid ) #ifdef SEAICE_CGRID CALL PLOT_FIELD_XYRL( seaiceMaskU, 'Current seaiceMaskU', & myIter, myThid ) CALL PLOT_FIELD_XYRL( seaiceMaskV, 'Current seaiceMaskV', & myIter, myThid ) #else CALL PLOT_FIELD_XYRL( UVM , 'Current UVM ' , & myIter, myThid ) #endif ENDIF #if (defined (SEAICE_CGRID) && defined (SEAICE_ALLOW_EVP)) DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx stressDivergenceX(I,J,bi,bj) = 0. _d 0 stressDivergenceY(I,J,bi,bj) = 0. _d 0 seaice_sigma1 (I,J,bi,bj) = 0. _d 0 seaice_sigma2 (I,J,bi,bj) = 0. _d 0 seaice_sigma12(I,J,bi,bj) = 0. _d 0 ENDDO ENDDO ENDDO ENDDO #endif /* SEAICE_ALLOW_EVP and SEAICE_CGRID */ C-- Set model variables to initial/restart conditions IF ( nIter0 .NE. 0 ) THEN CALL SEAICE_READ_PICKUP ( myThid ) ELSE DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx HSNOW(I,J,bi,bj)=0.2*HEFFM(i,j,bi,bj) YNEG(I,J,bi,bj)=ZERO TMIX(I,J,bi,bj)=TICE(I,J,bi,bj) DO k=1,3 HEFF(I,J,k,bi,bj)=SEAICE_initialHEFF*HEFFM(i,j,bi,bj) UICE(I,J,k,bi,bj)=ZERO VICE(I,J,k,bi,bj)=ZERO ENDDO ENDDO ENDDO ENDDO ENDDO C-- Read initial sea-ice thickness from file if available. IF ( HeffFile .NE. ' ' ) THEN _BEGIN_MASTER( myThid ) CALL READ_FLD_XY_RL( HeffFile, ' ', ZETA, 0, myThid ) _END_MASTER(myThid) _EXCH_XY_R8(ZETA,myThid) DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx DO k=1,3 HEFF(I,J,k,bi,bj) = MAX(ZETA(i,j,bi,bj),ZERO) ENDDO ENDDO ENDDO ENDDO ENDDO ENDIF DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx DO k=1,3 IF(HEFF(I,J,k,bi,bj).GT.ZERO) & AREA(I,J,k,bi,bj)=ONE ENDDO ENDDO ENDDO ENDDO ENDDO IF ( AreaFile .NE. ' ' ) THEN _BEGIN_MASTER( myThid ) CALL READ_FLD_XY_RL( AreaFile, ' ', ZETA, 0, myThid ) _END_MASTER(myThid) _EXCH_XY_R8(ZETA,myThid) DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx DO k=1,3 AREA(I,J,k,bi,bj) = MAX(ZETA(i,j,bi,bj),ZERO) AREA(I,J,k,bi,bj) = MIN(AREA(i,j,k,bi,bj),ONE) IF ( AREA(I,J,k,bi,bj) .LE. ZERO ) & HEFF(I,J,k,bi,bj) = ZERO IF ( HEFF(I,J,k,bi,bj) .LE. ZERO ) & AREA(I,J,k,bi,bj) = ZERO ENDDO ENDDO ENDDO ENDDO ENDDO ENDIF DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx IF(HEFF(I,J,1,bi,bj).LE.ZERO) & HSNOW(I,J,bi,bj)=ZERO ENDDO ENDDO ENDDO ENDDO ENDIF C--- Complete initialization DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx ZETA(I,J,bi,bj)=HEFF(I,J,1,bi,bj)*(1.0 _d 11) ETA(I,J,bi,bj)=ZETA(I,J,bi,bj)/4.0 _d 0 ENDDO ENDDO IF ( useRealFreshWaterFlux .AND. .NOT.useThSIce ) THEN DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx sIceLoad(i,j,bi,bj) = HEFF(I,J,1,bi,bj)*SEAICE_rhoIce & + HSNOW(I,J,bi,bj)* 330. _d 0 ENDDO ENDDO ENDIF ENDDO ENDDO RETURN END