C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/seaice/seaice_model.F,v 1.83 2010/11/18 17:32:37 jmc Exp $ C $Name: $ #include "SEAICE_OPTIONS.h" CBOP C !ROUTINE: SEAICE_MODEL C !INTERFACE: ========================================================== SUBROUTINE SEAICE_MODEL( myTime, myIter, myThid ) C !DESCRIPTION: \bv C *===========================================================* C | SUBROUTINE SEAICE_MODEL | C | o Time stepping of a dynamic/thermodynamic sea ice model. | C | Dynamics solver: Zhang/Hibler, JGR, 102, 8691-8702, 1997 | C | Thermodynamics: Hibler, MWR, 108, 1943-1973, 1980 | C | Rheology: Hibler, JPO, 9, 815- 846, 1979 | C | Snow: Zhang et al. , JPO, 28, 191- 217, 1998 | C | Parallel forward ice model written by Jinlun Zhang PSC/UW| C | & coupled into MITgcm by Dimitris Menemenlis (JPL) 2/2001| C | zhang@apl.washington.edu / menemenlis@jpl.nasa.gov | C *===========================================================* C *===========================================================* IMPLICIT NONE C \ev C !USES: =============================================================== #include "SIZE.h" #include "EEPARAMS.h" #include "DYNVARS.h" #include "PARAMS.h" #include "GRID.h" #include "FFIELDS.h" #include "SEAICE.h" #include "SEAICE_PARAMS.h" #ifdef ALLOW_EXF # include "EXF_OPTIONS.h" # include "EXF_FIELDS.h" #endif #ifdef ALLOW_SALT_PLUME # include "SALT_PLUME.h" #endif #ifdef ALLOW_AUTODIFF_TAMC # include "tamc.h" #endif C !INPUT PARAMETERS: =================================================== C myTime - Simulation time C myIter - Simulation timestep number C myThid - Thread no. that called this routine. _RL myTime INTEGER myIter INTEGER myThid CEndOfInterface C !LOCAL VARIABLES: ==================================================== C i,j,bi,bj :: Loop counters C iceFld :: Copy of seaice field INTEGER i, j, bi, bj c _RL iceFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy) CEOP #ifdef ALLOW_DEBUG IF (debugMode) CALL DEBUG_ENTER( 'SEAICE_MODEL', myThid ) #endif C-- Winds are from pkg/exf, which does not update edges. CALL EXCH_UV_AGRID_3D_RL( uwind, vwind, .TRUE., 1, myThid ) #ifdef ALLOW_THSICE IF ( useThSice ) THEN C-- Map thSice-variables to HEFF and AREA CALL SEAICE_MAP_THSICE( myTime, myIter, myThid ) ENDIF #endif /* ALLOW_THSICE */ IF ( .NOT.useThSice ) THEN #ifdef ALLOW_AUTODIFF_TAMC CADJ STORE heff = comlev1, key=ikey_dynamics, kind=isbyte CADJ STORE heffm = comlev1, key=ikey_dynamics, kind=isbyte CADJ STORE area = comlev1, key=ikey_dynamics, kind=isbyte CADJ STORE hsnow = comlev1, key=ikey_dynamics, kind=isbyte CADJ STORE tice = comlev1, key=ikey_dynamics, kind=isbyte #ifdef SEAICE_SALINITY CADJ STORE hsalt = comlev1, key=ikey_dynamics, kind=isbyte #endif #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,bi,bj).EQ.0.) & .OR.(area(i,j,bi,bj).EQ.0.) & ) THEN HEFF(i,j,bi,bj) = 0. _d 0 AREA(i,j,bi,bj) = 0. _d 0 HSNOW(i,j,bi,bj) = 0. _d 0 TICE(i,j,bi,bj) = celsius2K #ifdef SEAICE_SALINITY HSALT(i,j,bi,bj) = 0. _d 0 #endif ENDIF ENDDO ENDDO ENDDO ENDDO ENDIF #ifdef ALLOW_AUTODIFF_TAMC c DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=1-Oly,sNy+Oly DO i=1-Olx,sNx+Olx HEFFNM1(i,j,bi,bj) = 0. _d 0 AREANM1(i,j,bi,bj) = 0. _d 0 uIceNm1(i,j,bi,bj) = 0. _d 0 vIceNm1(i,j,bi,bj) = 0. _d 0 ENDDO ENDDO ENDDO ENDDO c CADJ STORE uwind = comlev1, key=ikey_dynamics, kind=isbyte CADJ STORE vwind = comlev1, key=ikey_dynamics, kind=isbyte CADJ STORE heff = comlev1, key=ikey_dynamics, kind=isbyte CADJ STORE heffm = comlev1, key=ikey_dynamics, kind=isbyte CADJ STORE area = comlev1, key=ikey_dynamics, kind=isbyte # ifdef SEAICE_ALLOW_DYNAMICS # ifdef SEAICE_CGRID CADJ STORE seaicemasku = comlev1, key=ikey_dynamics, kind=isbyte CADJ STORE seaicemaskv = comlev1, key=ikey_dynamics, kind=isbyte CADJ STORE fu = comlev1, key=ikey_dynamics, kind=isbyte CADJ STORE fv = comlev1, key=ikey_dynamics, kind=isbyte CADJ STORE uice = comlev1, key=ikey_dynamics, kind=isbyte CADJ STORE vice = comlev1, key=ikey_dynamics, kind=isbyte cphCADJ STORE eta = comlev1, key=ikey_dynamics, kind=isbyte cphCADJ STORE zeta = comlev1, key=ikey_dynamics, kind=isbyte cph( CADJ STORE dwatn = comlev1, key=ikey_dynamics, kind=isbyte cccCADJ STORE press0 = comlev1, key=ikey_dynamics, kind=isbyte cccCADJ STORE taux = comlev1, key=ikey_dynamics, kind=isbyte cccCADJ STORE tauy = comlev1, key=ikey_dynamics, kind=isbyte cccCADJ STORE zmax = comlev1, key=ikey_dynamics, kind=isbyte cccCADJ STORE zmin = comlev1, key=ikey_dynamics, kind=isbyte cph) # ifdef SEAICE_ALLOW_EVP CADJ STORE seaice_sigma1 = comlev1, key=ikey_dynamics, kind=isbyte CADJ STORE seaice_sigma2 = comlev1, key=ikey_dynamics, kind=isbyte CADJ STORE seaice_sigma12 = comlev1, key=ikey_dynamics, kind=isbyte # endif # endif # endif #endif /* ALLOW_AUTODIFF_TAMC */ C solve ice momentum equations and calculate ocean surface stress #ifdef ALLOW_DEBUG IF (debugMode) CALL DEBUG_CALL( 'SEAICE_DYNSOLVER', myThid ) #endif #ifdef SEAICE_CGRID CALL TIMER_START('SEAICE_DYNSOLVER [SEAICE_MODEL]',myThid) CALL SEAICE_DYNSOLVER ( myTime, myIter, myThid ) CALL TIMER_STOP ('SEAICE_DYNSOLVER [SEAICE_MODEL]',myThid) #else CALL TIMER_START('DYNSOLVER [SEAICE_MODEL]',myThid) CALL DYNSOLVER ( myTime, myIter, myThid ) CALL TIMER_STOP ('DYNSOLVER [SEAICE_MODEL]',myThid) #endif /* SEAICE_CGRID */ C-- Apply ice velocity open boundary conditions #ifdef ALLOW_OBCS # ifndef DISABLE_SEAICE_OBCS IF ( useOBCS ) CALL OBCS_APPLY_UVICE( uice, vice, myThid ) # endif /* DISABLE_SEAICE_OBCS */ #endif /* ALLOW_OBCS */ #ifdef ALLOW_THSICE IF ( .NOT.useThSice ) THEN #endif C-- Only call advection of heff, area, snow, and salt and C-- growth for the generic 0-layer thermodynamics of seaice C-- if useThSice=.false., otherwise the 3-layer Winton thermodynamics C-- (called from DO_OCEANIC_PHYSICS) take care of this C NOW DO ADVECTION and DIFFUSION IF ( SEAICEadvHeff .OR. SEAICEadvArea .OR. SEAICEadvSnow & .OR. SEAICEadvSalt .OR. SEAICEadvAge ) THEN #ifdef ALLOW_DEBUG IF (debugMode) CALL DEBUG_CALL( 'SEAICE_ADVDIFF', myThid ) #endif CALL SEAICE_ADVDIFF( myTime, myIter, myThid ) ENDIF #ifdef ALLOW_AUTODIFF_TAMC CADJ STORE heffm = comlev1, key=ikey_dynamics, kind=isbyte cph-test( cphCADJ STORE heff = comlev1, key=ikey_dynamics, kind=isbyte cphCADJ STORE area = comlev1, key=ikey_dynamics, kind=isbyte cphCADJ STORE hsnow = comlev1, key=ikey_dynamics, kind=isbyte cphCADJ STORE qnet = comlev1, key=ikey_dynamics, kind=isbyte cphCADJ STORE qsw = comlev1, key=ikey_dynamics, kind=isbyte cphCADJ STORE tice = comlev1, key=ikey_dynamics, kind=isbyte cph-test) # ifdef SEAICE_ALLOW_DYNAMICS cphCADJ STORE uice = comlev1, key=ikey_dynamics, kind=isbyte cphCADJ STORE vice = comlev1, key=ikey_dynamics, kind=isbyte # endif #endif /* ALLOW_AUTODIFF_TAMC */ C thermodynamics growth C must call growth after calling advection C because of ugly time level business IF ( usePW79thermodynamics ) THEN #ifdef ALLOW_DEBUG IF (debugMode) CALL DEBUG_CALL( 'SEAICE_GROWTH', myThid ) #endif #ifndef SEAICE_ALLOW_TD_IF CALL SEAICE_GROWTH( myTime, myIter, myThid ) #else CALL SEAICE_GROWTH_IF( myTime, myIter, myThid ) #endif ENDIF C-- Apply ice tracer open boundary conditions #ifdef ALLOW_OBCS # ifndef DISABLE_SEAICE_OBCS IF ( useOBCS ) CALL OBCS_APPLY_SEAICE( myThid ) # endif /* DISABLE_SEAICE_OBCS */ #endif /* ALLOW_OBCS */ C-- Update overlap regions for a bunch of stuff _EXCH_XY_RL( HEFF, myThid ) _EXCH_XY_RL( AREA, myThid ) _EXCH_XY_RL( HSNOW, myThid ) #ifdef SEAICE_SALINITY _EXCH_XY_RL( HSALT, myThid ) #endif #ifdef SEAICE_AGE _EXCH_XY_RL( IceAge,myThid ) #endif _EXCH_XY_RS(EmPmR, myThid ) _EXCH_XY_RS(saltFlux, myThid ) _EXCH_XY_RS(Qnet , myThid ) #ifdef SHORTWAVE_HEATING _EXCH_XY_RS(Qsw , myThid ) #endif /* SHORTWAVE_HEATING */ #ifdef ALLOW_SALT_PLUME IF ( useSALT_PLUME ) & _EXCH_XY_RL(saltPlumeFlux, myThid ) #endif /* ALLOW_SALT_PLUME */ #ifdef ATMOSPHERIC_LOADING IF ( useRealFreshWaterFlux ) & _EXCH_XY_RS( sIceLoad, myThid ) #endif C-- In case we use scheme with a large stencil that extends into overlap: #ifdef ALLOW_OBCS IF ( useOBCS ) THEN DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) CALL OBCS_COPY_TRACER( HEFF(1-Olx,1-Oly,bi,bj), I 1, bi, bj, myThid ) CALL OBCS_COPY_TRACER( AREA(1-Olx,1-Oly,bi,bj), I 1, bi, bj, myThid ) CALL OBCS_COPY_TRACER( HSNOW(1-Olx,1-Oly,bi,bj), I 1, bi, bj, myThid ) #ifdef SEAICE_SALINITY CALL OBCS_COPY_TRACER( HSALT(1-Olx,1-Oly,bi,bj), I 1, bi, bj, myThid ) #endif #ifdef SEAICE_AGE CALL OBCS_COPY_TRACER( IceAge(1-Olx,1-Oly,bi,bj), I 1, bi, bj, myThid ) #endif ENDDO ENDDO ENDIF #endif /* ALLOW_OBCS */ #ifdef ALLOW_DIAGNOSTICS IF ( useDiagnostics ) THEN C diagnostics for "non-state variables" that are modified by C the seaice model # ifdef ALLOW_EXF CALL DIAGNOSTICS_FILL(UWIND ,'SIuwind ',0,1 ,0,1,1,myThid) CALL DIAGNOSTICS_FILL(VWIND ,'SIvwind ',0,1 ,0,1,1,myThid) # endif CALL DIAGNOSTICS_FILL_RS(FU ,'SIfu ',0,1 ,0,1,1,myThid) CALL DIAGNOSTICS_FILL_RS(FV ,'SIfv ',0,1 ,0,1,1,myThid) CALL DIAGNOSTICS_FILL_RS(EmPmR,'SIempmr ',0,1 ,0,1,1,myThid) CALL DIAGNOSTICS_FILL_RS(Qnet ,'SIqnet ',0,1 ,0,1,1,myThid) CALL DIAGNOSTICS_FILL_RS(Qsw ,'SIqsw ',0,1 ,0,1,1,myThid) ENDIF #endif /* ALLOW_DIAGNOSTICS */ #ifdef ALLOW_THSICE C endif .not.useThSice ENDIF #endif /* ALLOW_THSICE */ CML This has already been done in seaice_ocean_stress/ostres, so why repeat? CML CALL EXCH_UV_XY_RS(fu,fv,.TRUE.,myThid) #ifdef ALLOW_EXF # ifdef ALLOW_AUTODIFF_TAMC # if (defined (ALLOW_AUTODIFF_MONITOR)) CALL EXF_ADJOINT_SNAPSHOTS( 3, myTime, myIter, myThid ) # endif # endif #endif #ifdef ALLOW_DEBUG IF (debugMode) CALL DEBUG_LEAVE( 'SEAICE_MODEL', myThid ) #endif RETURN END