pkg/thsice/thsice_main.F

fc7306ba7d Jean*0001 #include "THSICE_OPTIONS.h"
6b47d550f4 Mart*0002 #ifdef ALLOW_AUTODIFF
                0003 # include "AUTODIFF_OPTIONS.h"
c1c3d0f9d7 Patr*0004 # ifdef ALLOW_EXF
                0005 #  include "EXF_OPTIONS.h"
                0006 # endif
                0007 #endif
7269783f6f Jean*0008
87ea84cac6 Jean*0009 CBOP
fc7306ba7d Jean*0010 C     !ROUTINE: THSICE_MAIN
                0011 C     !INTERFACE:
7269783f6f Jean*0012       SUBROUTINE THSICE_MAIN(
fc7306ba7d Jean*0013      I                        myTime, myIter, myThid )
87ea84cac6 Jean*0014 C     !DESCRIPTION: \bv
fc7306ba7d Jean*0015 C     *==========================================================*
7269783f6f Jean*0016 C     | S/R  THSICE_MAIN
                0017 C     | o Therm_SeaIce main routine.
fc7306ba7d Jean*0018 C     |   step forward Thermodynamic_SeaIce variables and modify
                0019 C     |    ocean surface forcing accordingly.
                0020 C     *==========================================================*
                0021
                0022 C     !USES:
                0023       IMPLICIT NONE
87ea84cac6 Jean*0024
fc7306ba7d Jean*0025 C     === Global variables ===
                0026 #include "SIZE.h"
                0027 #include "EEPARAMS.h"
                0028 #include "PARAMS.h"
                0029 #include "FFIELDS.h"
                0030 #include "THSICE_PARAMS.h"
2c032d7179 Gael*0031 #include "THSICE_SIZE.h"
87ea84cac6 Jean*0032 #include "THSICE_VARS.h"
6b47d550f4 Mart*0033 #ifdef ALLOW_AUTODIFF
9fbb8d18a8 Jean*0034 # include "THSICE_COST.h"
9439f3829d Jean*0035 # include "DYNVARS.h"
c1c3d0f9d7 Patr*0036 # ifdef ALLOW_EXF
                0037 #  include "EXF_PARAM.h"
                0038 #  include "EXF_CONSTANTS.h"
6b47d550f4 Mart*0039 #  include "EXF_FIELDS.h"
c1c3d0f9d7 Patr*0040 # endif /* ALLOW_EXF */
6b47d550f4 Mart*0041 #endif /* ALLOW_AUTODIFF */
                0042 #ifdef ALLOW_AUTODIFF_TAMC
                0043 # include "tamc.h"
c567874792 Patr*0044 #endif
7269783f6f Jean*0045
fc7306ba7d Jean*0046 C     !INPUT/OUTPUT PARAMETERS:
                0047 C     === Routine arguments ===
a4eca6e929 Jean*0048 C     myTime    :: Current time in simulation (s)
                0049 C     myIter    :: Current iteration number
                0050 C     myThid    :: My Thread Id. number
                0051       _RL     myTime
fc7306ba7d Jean*0052       INTEGER myIter
                0053       INTEGER myThid
87ea84cac6 Jean*0054 CEOP
fc7306ba7d Jean*0055
                0056 #ifdef ALLOW_THSICE
                0057 C     !LOCAL VARIABLES:
                0058 C     === Local variables ===
b947c8cb90 Jean*0059 C     prcAtm    :: total precip from the atmosphere [kg/m2/s]
                0060 C     snowPr    :: snow precipitation               [kg/m2/s]
bd7be113e1 Jean*0061 C     qPrcRn    :: Energy content of Precip+RunOff (+=down) [W/m2]
fc7306ba7d Jean*0062       INTEGER i,j
                0063       INTEGER bi,bj
                0064       INTEGER iMin, iMax
                0065       INTEGER jMin, jMax
de836be2bc Jean*0066       _RL prcAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
c8bafbe9ad Jean*0067       _RL snowPr(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
bd7be113e1 Jean*0068       _RL qPrcRn(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
7269783f6f Jean*0069 c     _RL evpAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0070 c     _RL flxAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0071 c     _RL flxSW (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
fc7306ba7d Jean*0072       _RL tauFac
40d541aac0 Jean*0073 #ifdef ALLOW_EXF
                0074       INTEGER grpDiag
                0075 #endif
7c50f07931 Mart*0076 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0077 C     tkey  :: tape key (depends on tiles)
                0078       INTEGER tkey
7c50f07931 Mart*0079 #endif
fc7306ba7d Jean*0080
                0081 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0082
40d541aac0 Jean*0083 #ifdef ALLOW_EXF
                0084       IF ( useEXF .AND. useDiagnostics ) THEN
                0085 C-    Fill-in EXF surface flux diags, weighted by open-ocean fraction
                0086         grpDiag = 2
                0087         IF ( thSIce_skipThermo ) grpDiag = -2
                0088         CALL EXF_WEIGHT_SFX_DIAGS(
                0089      I                  iceMask, grpDiag, myTime, myIter, myThid )
a3aa8e4116 Jean*0090         IF ( .NOT.useSEAICE ) CALL EXF_WEIGHT_SFX_DIAGS(
                0091      I                       iceMask, -1, myTime, myIter, myThid )
40d541aac0 Jean*0092       ENDIF
                0093 #endif /* ALLOW_EXF */
                0094
24d0c27edd Jean*0095 C-    Only compute/update seaice fields over the interior
                0096 C     (excluding overlap) and apply exchanges when needed
                0097       iMin = 1
                0098       iMax = sNx
                0099       jMin = 1
                0100       jMax = sNy
fc7306ba7d Jean*0101
                0102       DO bj=myByLo(myThid),myByHi(myThid)
                0103        DO bi=myBxLo(myThid),myBxHi(myThid)
87ea84cac6 Jean*0104
c567874792 Patr*0105 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0106         tkey = bi + (bj-1)*nSx + (ikey_dynamics-1)*nSx*nSy
                0107 CADJ STORE uvel  (:,:,1,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
                0108 CADJ STORE vvel  (:,:,1,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
                0109 CADJ STORE qsw     (:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
                0110 CADJ STORE Qice1   (:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
                0111 CADJ STORE Qice2   (:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
                0112 CADJ STORE sHeating(:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
                0113 CADJ STORE Tsrf    (:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
c567874792 Patr*0114 #endif
                0115
9439f3829d Jean*0116         DO j=1-OLy,sNy+OLy
                0117          DO i=1-OLx,sNx+OLx
de836be2bc Jean*0118           prcAtm  (i,j,bi,bj) = 0. _d 0
c8bafbe9ad Jean*0119           snowPr  (i,j) = 0. _d 0
bd7be113e1 Jean*0120           qPrcRn  (i,j) = 0. _d 0
87ea84cac6 Jean*0121          ENDDO
6b6ed88e13 Mart*0122         ENDDO
c567874792 Patr*0123
6b47d550f4 Mart*0124 #ifndef ALLOW_AUTODIFF
9fbb8d18a8 Jean*0125         IF ( .NOT.useCheapAML ) THEN
                0126 #endif
147f2f3fa7 Jean*0127          CALL THSICE_GET_OCEAN(
                0128      I                          bi, bj, myTime, myIter, myThid )
9439f3829d Jean*0129
c567874792 Patr*0130 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0131 CADJ STORE iceMask   (:,:,bi,bj) = comlev1_bibj,key=tkey,kind=isbyte
                0132 CADJ STORE iceHeight (:,:,bi,bj) = comlev1_bibj,key=tkey,kind=isbyte
                0133 CADJ STORE snowHeight(:,:,bi,bj) = comlev1_bibj,key=tkey,kind=isbyte
                0134 CADJ STORE snowAge   (:,:,bi,bj) = comlev1_bibj,key=tkey,kind=isbyte
                0135 CADJ STORE hOceMxL   (:,:,bi,bj) = comlev1_bibj,key=tkey,kind=isbyte
                0136 CADJ STORE tOceMxL   (:,:,bi,bj) = comlev1_bibj,key=tkey,kind=isbyte
                0137 CADJ STORE sOceMxL   (:,:,bi,bj) = comlev1_bibj,key=tkey,kind=isbyte
                0138 CADJ STORE v2ocMxL   (:,:,bi,bj) = comlev1_bibj,key=tkey,kind=isbyte
6b47d550f4 Mart*0139 #endif /* ALLOW_AUTODIFF_TAMC */
                0140
                0141 #ifndef ALLOW_AUTODIFF
147f2f3fa7 Jean*0142 C-   end if not useCheapAML
                0143         ENDIF
c567874792 Patr*0144 #endif
                0145
6b47d550f4 Mart*0146 #ifndef ALLOW_AUTODIFF
147f2f3fa7 Jean*0147         IF ( useBulkforce .OR. useCheapAML ) THEN
6b6ed88e13 Mart*0148          CALL THSICE_GET_PRECIP(
e22ef29cc8 Jean*0149      I                  iceMask, tOceMxL,
de836be2bc Jean*0150      O                  prcAtm(1-OLx,1-OLy,bi,bj),
c8bafbe9ad Jean*0151      O                  snowPr, qPrcRn,
7269783f6f Jean*0152      O                  icFlxSW(1-OLx,1-OLy,bi,bj),
2d104ee1fd Jean*0153      I                  iMin,iMax,jMin,jMax, bi,bj, myThid )
6b6ed88e13 Mart*0154         ENDIF
2d104ee1fd Jean*0155 #endif
6b6ed88e13 Mart*0156         IF ( useEXF ) THEN
                0157          CALL THSICE_MAP_EXF(
bd7be113e1 Jean*0158      I                  iceMask, tOceMxL,
de836be2bc Jean*0159      O                  prcAtm(1-OLx,1-OLy,bi,bj),
c8bafbe9ad Jean*0160      O                  snowPr, qPrcRn,
2a9474d935 Mart*0161      O                  icFlxSW(1-OLx,1-OLy,bi,bj),
                0162      I                  iMin,iMax,jMin,jMax, bi,bj, myThid )
6b6ed88e13 Mart*0163         ENDIF
87ea84cac6 Jean*0164
6b47d550f4 Mart*0165 #ifndef ALLOW_AUTODIFF
147f2f3fa7 Jean*0166         IF ( .NOT.( useCheapAML .OR. thSIce_skipThermo ) ) THEN
1818702d6f Patr*0167 #endif
147f2f3fa7 Jean*0168          CALL THSICE_STEP_TEMP(
7269783f6f Jean*0169      I                     bi, bj, iMin, iMax, jMin, jMax,
                0170      I                     myTime, myIter, myThid )
                0171
d6f06800ae Patr*0172 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0173 C     These store directives avoid recomputation of this routine in ad-code.
                0174 CADJ STORE flxCndBt(:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
                0175 CADJ STORE EmPmR   (:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
                0176 CADJ STORE Qnet    (:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
                0177 CADJ STORE snowAge (:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
                0178 CADJ STORE sHeating(:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
                0179 CADJ STORE snowPr              = comlev1_bibj, key=tkey, kind=isbyte
                0180 CADJ STORE Qice1   (:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
                0181 CADJ STORE Qice2   (:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
                0182 CADJ STORE icFrwAtm(:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
a85293d087 Mart*0183 C     This part of the state is generally not necessary, but you never
edb6656069 Mart*0184 C     know so I leave it, but commented out.
                0185 cCADJ STORE Qsw     (:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
                0186 cCADJ STORE icFlxSW (:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
                0187 cCADJ STORE Tice1   (:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
                0188 cCADJ STORE Tice2   (:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
                0189 cCADJ STORE sHeating(:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
                0190 cCADJ STORE Tsrf    (:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
                0191 cCADJ STORE icFlxAtm(:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
6b47d550f4 Mart*0192 #endif
                0193
                0194 #ifndef ALLOW_AUTODIFF
147f2f3fa7 Jean*0195 C-   end if not skipThermo / useCheapAML
9439f3829d Jean*0196         ENDIF
                0197         IF ( .NOT.thSIce_skipThermo ) THEN
d6f06800ae Patr*0198 #endif
147f2f3fa7 Jean*0199          CALL THSICE_STEP_FWD(
7269783f6f Jean*0200      I                     bi, bj, iMin, iMax, jMin, jMax,
c8bafbe9ad Jean*0201      I                     prcAtm(1-OLx,1-OLy,bi,bj),
                0202      I                     snowPr, qPrcRn,
fc7306ba7d Jean*0203      I                     myTime, myIter, myThid )
6b47d550f4 Mart*0204 #ifndef ALLOW_AUTODIFF
3772953427 Jean*0205         ELSE
                0206 C-    Compute sIceLoad (if not previously done)
                0207           DO j=1,sNy
                0208            DO i=1,sNx
                0209              sIceLoad(i,j,bi,bj) = ( snowHeight(i,j,bi,bj)*rhos
                0210      &                             + iceHeight(i,j,bi,bj)*rhoi
                0211      &                             )*iceMask(i,j,bi,bj)
                0212            ENDDO
                0213           ENDDO
8a23e1b5d8 Jean*0214         ENDIF
                0215 #endif
87ea84cac6 Jean*0216
de836be2bc Jean*0217 C--  end bi,bj loop
                0218        ENDDO
                0219       ENDDO
                0220
                0221 #ifdef ALLOW_BALANCE_FLUXES
                0222 C--   Balance net Fresh-Water flux from Atm+Land
                0223       IF ( thSIceBalanceAtmFW.NE.0 ) THEN
8a23e1b5d8 Jean*0224         CALL THSICE_BALANCE_FRW(
de836be2bc Jean*0225      I                      iMin, iMax, jMin, jMax,
                0226      I                      prcAtm, myTime, myIter, myThid )
                0227       ENDIF
6b47d550f4 Mart*0228 #endif /* ALLOW_BALANCE_FLUXES */
de836be2bc Jean*0229
e3dfb30901 Jean*0230 C     add a small piece of code to check AddFluid implementation:
                0231 c#include "thsice_test_addfluid.h"
                0232
                0233 C--   Exchange fields that are advected by seaice dynamics
96502f6244 Jean*0234       IF ( useSEAICE .OR. thSIceAdvScheme.GT.0
24d0c27edd Jean*0235      &               .OR. stressReduction.GT.zeroRL ) THEN
96502f6244 Jean*0236         CALL THSICE_DO_EXCH( myThid )
e3dfb30901 Jean*0237       ENDIF
3772953427 Jean*0238       IF ( thSIceAdvScheme.GT.0 .AND. .NOT.useSEAICE ) THEN
96502f6244 Jean*0239 C-    when useSEAICE=.true., this S/R is called from SEAICE_MODEL;
                0240 C     otherwise, call it from here, after thsice-thermodynamics is done
                0241          CALL THSICE_DO_ADVECT(
                0242      I                          0, 0, myTime, myIter, myThid )
3772953427 Jean*0243       ELSEIF ( thSIceAdvScheme.LE.0 .AND. useRealFreshWaterFlux ) THEN
                0244         _EXCH_XY_RS( sIceLoad, myThid )
96502f6244 Jean*0245       ENDIF
e3dfb30901 Jean*0246
                0247       DO bj=myByLo(myThid),myByHi(myThid)
                0248        DO bi=myBxLo(myThid),myBxHi(myThid)
2f1d6bb332 Patr*0249 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0250         tkey = bi + (bj-1)*nSx + (ikey_dynamics-1)*nSx*nSy
                0251 CADJ STORE iceMask(:,:,bi,bj) = comlev1_bibj,key=tkey,kind=isbyte
                0252 CADJ STORE fu     (:,:,bi,bj) = comlev1_bibj,key=tkey,kind=isbyte
                0253 CADJ STORE fv     (:,:,bi,bj) = comlev1_bibj,key=tkey,kind=isbyte
2f1d6bb332 Patr*0254 #endif /* ALLOW_AUTODIFF_TAMC */
                0255
c9573f2063 Jean*0256 C--   Cumulate time-averaged fields and also fill-up flux diagnostics
                0257 C     (if not done in THSICE_DO_ADVECT call)
                0258         IF ( thSIceAdvScheme.LE.0 ) THEN
                0259          CALL THSICE_AVE(
                0260      I                     bi,bj, myTime, myIter, myThid )
                0261         ENDIF
b7ebc7bc01 Jean*0262 C--   note: If useSEAICE=.true., the stress is computed in seaice_model,
                0263 C--   and stressReduction is always set to zero
                0264         IF ( stressReduction.GT. 0. _d 0 ) THEN
24d0c27edd Jean*0265           DO j = 1-OLy,sNy+OLy-1
                0266            DO i = 2-OLx,sNx+OLx-1
fc7306ba7d Jean*0267             tauFac = stressReduction
                0268      &             *(iceMask(i-1,j,bi,bj)+iceMask(i,j,bi,bj))*0.5 _d 0
                0269             fu(i,j,bi,bj) = (1. _d 0 - tauFac)*fu(i,j,bi,bj)
6b6ed88e13 Mart*0270            ENDDO
fc7306ba7d Jean*0271           ENDDO
24d0c27edd Jean*0272           DO j = 2-OLy,sNy+OLy-1
                0273            DO i = 1-OLx,sNx+OLx-1
fc7306ba7d Jean*0274             tauFac = stressReduction
                0275      &             *(iceMask(i,j-1,bi,bj)+iceMask(i,j,bi,bj))*0.5 _d 0
                0276             fv(i,j,bi,bj) = (1. _d 0 - tauFac)*fv(i,j,bi,bj)
6b6ed88e13 Mart*0277            ENDDO
fc7306ba7d Jean*0278           ENDDO
                0279         ENDIF
                0280
                0281 C--  end bi,bj loop
                0282        ENDDO
                0283       ENDDO
                0284
                0285 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
6b47d550f4 Mart*0286 #endif /* ALLOW_THSICE */
fc7306ba7d Jean*0287
                0288       RETURN
                0289       END