pkg/thsice/thsice_step_temp.F

61c3f3d76b Jean*0001 #include "THSICE_OPTIONS.h"
6b47d550f4 Mart*0002 #ifdef ALLOW_AUTODIFF
                0003 # include "AUTODIFF_OPTIONS.h"
                0004 #endif
61c3f3d76b Jean*0005
                0006 CBOP
                0007 C     !ROUTINE: THSICE_STEP_TEMP
                0008 C     !INTERFACE:
                0009       SUBROUTINE THSICE_STEP_TEMP(
                0010      I             bi, bj, iMin, iMax, jMin, jMax,
                0011      I             myTime, myIter, myThid )
                0012 C     !DESCRIPTION: \bv
                0013 C     *==========================================================*
                0014 C     | S/R  THSICE_STEP_TEMP
                0015 C     | o Step Forward Surface and SeaIce Temperature.
                0016 C     *==========================================================*
                0017 C     \ev
                0018
                0019 C     !USES:
                0020       IMPLICIT NONE
                0021
                0022 C     === Global variables ===
                0023 #include "SIZE.h"
                0024 #include "EEPARAMS.h"
                0025 #include "PARAMS.h"
                0026 #include "FFIELDS.h"
                0027 #include "THSICE_SIZE.h"
                0028 #include "THSICE_PARAMS.h"
                0029 #include "THSICE_VARS.h"
d6f06800ae Patr*0030 #ifdef ALLOW_AUTODIFF_TAMC
                0031 # include "tamc.h"
                0032 #endif
                0033
61c3f3d76b Jean*0034       INTEGER siLo, siHi, sjLo, sjHi
                0035       PARAMETER ( siLo = 1-OLx , siHi = sNx+OLx )
                0036       PARAMETER ( sjLo = 1-OLy , sjHi = sNy+OLy )
                0037
                0038 C     !INPUT/OUTPUT PARAMETERS:
                0039 C     === Routine arguments ===
                0040 C- input:
                0041 C     bi,bj   :: tile indices
                0042 C   iMin,iMax :: computation domain: 1rst index range
                0043 C   jMin,jMax :: computation domain: 2nd  index range
0b1ccf0764 Jean*0044 C     myTime  :: Current time in simulation
                0045 C     myIter  :: Current iteration number
                0046 C     myThid  :: my Thread Id number
61c3f3d76b Jean*0047 C-- Modify fluxes hold in commom blocks
                0048 C- input:
                0049 C     icFlxSW :: (Inp) short-wave heat flux (+=down): downward comp. only
                0050 C- output
                0051 C     icFlxSW :: (Out) net SW flux into ocean (+=down)
                0052 C     icFlxAtm:: net flux of energy from the atmosphere [W/m2] (+=down)
                0053 C     icFrwAtm:: evaporation to the atmosphere (kg/m2/s) (>0 if evaporate)
                0054 C--
                0055       INTEGER bi,bj
                0056       INTEGER iMin, iMax
                0057       INTEGER jMin, jMax
                0058       _RL  myTime
                0059       INTEGER myIter
                0060       INTEGER myThid
                0061 CEOP
                0062
                0063 #ifdef ALLOW_THSICE
                0064 C     !LOCAL VARIABLES:
                0065 C     === Local variables ===
                0066 C     tFrzOce   :: sea-water freezing temperature [oC] (function of S)
                0067 C     dTsrf     :: surf. temp adjusment: Ts^n+1 - Ts^n
b022b1e505 Jean*0068 C     tmpFlx    :: dummy array for surface fluxes and derivative vs Tsurf
                0069 C Note: dTsrf & tmpFlx are not used here; just allocate enough space for dTsrf.
61c3f3d76b Jean*0070       INTEGER i,j
                0071       _RL tFrzOce(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
a85293d087 Mart*0072       _RL dTsrf  (1:sNx,1:sNy)
                0073 C     some compilers complain about the wrong dimensions of this field
                0074 C     so it has to be defined with overlaps
                0075       _RL tmpFlx (1:sNx,1:sNy,0:2)
61c3f3d76b Jean*0076       _RL opFrac, icFrac
                0077       LOGICAL dBugFlag
                0078
7c50f07931 Mart*0079 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0080 C     tkey  :: tape key (depends on tiles)
                0081       INTEGER tkey
7c50f07931 Mart*0082 #endif
61c3f3d76b Jean*0083 C-    define grid-point location where to print debugging values
                0084 #include "THSICE_DEBUG.h"
                0085
                0086  1010 FORMAT(A,1P4E14.6)
                0087
                0088 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0089
d6f06800ae Patr*0090 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0091       tkey = bi + (bj-1)*nSx + (ikey_dynamics-1)*nSx*nSy
6b47d550f4 Mart*0092 #endif /* ALLOW_AUTODIFF_TAMC */
                0093
a85293d087 Mart*0094 C     Initialise some local arrays
c1c3d0f9d7 Patr*0095       DO j = 1-OLy, sNy+OLy
                0096        DO i = 1-OLx, sNx+OLx
a85293d087 Mart*0097         tFrzOce(i,j)  = 0. _d 0
                0098        ENDDO
                0099       ENDDO
                0100       DO j = 1, sNy
                0101        DO i = 1, sNx
                0102         dTsrf(i,j)    = 0. _d 0
                0103         tmpFlx(i,j,0) = 0. _d 0
                0104         tmpFlx(i,j,1) = 0. _d 0
                0105         tmpFlx(i,j,2) = 0. _d 0
c1c3d0f9d7 Patr*0106        ENDDO
                0107       ENDDO
d6f06800ae Patr*0108
ae605e558b Jean*0109       dBugFlag = debugLevel.GE.debLevC
61c3f3d76b Jean*0110 C-    Initialise flxAtm,evpAtm
a85293d087 Mart*0111       DO j = 1-OLy, sNy+OLy
                0112        DO i = 1-OLx, sNx+OLx
                0113         icFlxAtm(i,j,bi,bj) = 0.
                0114         icFrwAtm(i,j,bi,bj) = 0.
61c3f3d76b Jean*0115        ENDDO
a85293d087 Mart*0116       ENDDO
61c3f3d76b Jean*0117
d6f06800ae Patr*0118 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0119 CADJ STORE tsrf(:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
d6f06800ae Patr*0120 #endif
61c3f3d76b Jean*0121 c     IF ( fluidIsWater ) THEN
                0122        CALL THSICE_ALBEDO(
                0123      I          bi, bj, siLo, siHi, sjLo, sjHi,
                0124      I          iMin,iMax, jMin,jMax,
                0125      I          iceMask(siLo,sjLo,bi,bj), iceHeight(siLo,sjLo,bi,bj),
                0126      I          snowHeight(siLo,sjLo,bi,bj), Tsrf(siLo,sjLo,bi,bj),
                0127      I          snowAge(siLo,sjLo,bi,bj),
ce354ad541 Jean*0128      O          siceAlb(siLo,sjLo,bi,bj), icAlbNIR(siLo,sjLo,bi,bj),
61c3f3d76b Jean*0129      I          myTime, myIter, myThid )
                0130
                0131 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0132 C    part.1 : ice-covered fraction ;
                0133 C     Solve for surface and ice temperature (implicitly) ; compute surf. fluxes
                0134 C-------
d6f06800ae Patr*0135
                0136 #ifdef ALLOW_AUTODIFF_TAMC
a85293d087 Mart*0137 C     avoid calling s/r thsice_albedosolve4temp in adjoint routine again
edb6656069 Mart*0138 CADJ STORE siceAlb (:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
                0139 CADJ STORE icAlbNIR(:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
                0140 CADJ STORE icFlxSW (:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
d6f06800ae Patr*0141 #endif
                0142
9a36a81763 Mart*0143 #ifdef ALLOW_DBUG_THSICE
61c3f3d76b Jean*0144        DO j = jMin, jMax
                0145         DO i = iMin, iMax
                0146          IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN
                0147           IF ( dBug(i,j,bi,bj) ) THEN
                0148            WRITE(6,'(A,2I4,2I2)') 'ThSI_STEP_T: i,j=',i,j,bi,bj
                0149            WRITE(6,1010) 'ThSI_STEP_T: iceMask, hIc, hSn, Tsf  =',
                0150      &                   iceMask(i,j,bi,bj), iceHeight(i,j,bi,bj),
                0151      &                   snowHeight(i,j,bi,bj), Tsrf(i,j,bi,bj)
                0152           ENDIF
9a36a81763 Mart*0153          ENDIF
                0154         ENDDO
                0155        ENDDO
61c3f3d76b Jean*0156 #endif
9a36a81763 Mart*0157        DO j = jMin, jMax
                0158         DO i = iMin, iMax
                0159          IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN
61c3f3d76b Jean*0160 C-      surface net SW flux:
                0161           icFlxSW(i,j,bi,bj) = icFlxSW(i,j,bi,bj)
                0162      &                       *(1. _d 0 - siceAlb(i,j,bi,bj))
                0163           tFrzOce(i,j) = -mu_Tf*sOceMxL(i,j,bi,bj)
                0164          ELSE
                0165           tFrzOce(i,j) = 0. _d 0
                0166          ENDIF
                0167         ENDDO
                0168        ENDDO
                0169
77008a74ba Patr*0170 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0171 CADJ STORE Qice1   (:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
                0172 CADJ STORE Qice2   (:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
                0173 CADJ STORE sHeating(:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
                0174 CADJ STORE Tsrf    (:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
77008a74ba Patr*0175 #endif
61c3f3d76b Jean*0176        CALL THSICE_SOLVE4TEMP(
6dc8890c80 Patr*0177      I          bi, bj,
b022b1e505 Jean*0178      I          iMin,iMax, jMin,jMax, dBugFlag,
2a9474d935 Mart*0179      I          useBulkForce, useEXF,
61c3f3d76b Jean*0180      I          iceMask(siLo,sjLo,bi,bj), iceHeight(siLo,sjLo,bi,bj),
b022b1e505 Jean*0181      I          snowHeight(siLo,sjLo,bi,bj), tFrzOce, tmpFlx,
61c3f3d76b Jean*0182      U          icFlxSW(siLo,sjLo,bi,bj), Tsrf(siLo,sjLo,bi,bj),
                0183      U          Qice1(siLo,sjLo,bi,bj), Qice2(siLo,sjLo,bi,bj),
b022b1e505 Jean*0184      O          Tice1(siLo,sjLo,bi,bj), Tice2(siLo,sjLo,bi,bj), dTsrf,
61c3f3d76b Jean*0185      O          sHeating(siLo,sjLo,bi,bj), flxCndBt(siLo,sjLo,bi,bj),
                0186      O          icFlxAtm(siLo,sjLo,bi,bj), icFrwAtm(siLo,sjLo,bi,bj),
                0187      I          myTime, myIter, myThid )
                0188
a85293d087 Mart*0189 #ifdef ALLOW_AUTODIFF_TAMC
                0190 C     avoid calling s/r thsice_solve4temp in adjoint routine again
edb6656069 Mart*0191 CADJ STORE icFlxSW(:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
a85293d087 Mart*0192 #endif
61c3f3d76b Jean*0193        DO j = jMin, jMax
                0194         DO i = iMin, iMax
                0195          IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN
                0196           icFrac  = iceMask(i,j,bi,bj)
                0197           opFrac = 1. _d 0 - icFrac
                0198 C--    Update Fluxes :
                0199           Qsw(i,j,bi,bj) = opFrac*Qsw(i,j,bi,bj)
                0200      &                   - icFrac*icFlxSW(i,j,bi,bj)
                0201          ENDIF
                0202         ENDDO
                0203        ENDDO
                0204 c     ENDIF
                0205
                0206 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0207 #endif /* ALLOW_THSICE */
                0208
                0209       RETURN
                0210       END