Back to home page

MITgcm

 
 

    


File indexing completed on 2018-03-02 18:38:25 UTC

view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
8fd83faf35 Jean*0001 #include "CHEAPAML_OPTIONS.h"
                0002 #ifdef ALLOW_THSICE
                0003 # include "THSICE_OPTIONS.h"
                0004 #endif
                0005 #ifdef ALLOW_SEAICE
                0006 # include "SEAICE_OPTIONS.h"
                0007 #endif
                0008 
                0009 CBOP
                0010 C     !ROUTINE: CHEAPAML_SEAICE
                0011 C     !INTERFACE:
                0012       SUBROUTINE CHEAPAML_SEAICE(
                0013      I                    swDown, lwDown, uWind, vWind, LVapor,
                0014      O                    fsha, flha, evp, xolw, ssqt, q100, cdq,
                0015      O                    Tsurf, iceFrac, sw2oce,
                0016      I                    bi, bj, myTime, myIter, myThid )
                0017 C     !DESCRIPTION: \bv
                0018 C     *==========================================================*
                0019 C     | S/R CHEAPAML_SEAICE
                0020 C     | o Compute fluxes over seaice by calling seaice routine
                0021 C     |   to solve for surface temperature.
                0022 C     *==========================================================*
                0023 C     \ev
                0024 
                0025 C     !USES:
                0026       IMPLICIT NONE
                0027 C     == Global variables ===
                0028 #include "SIZE.h"
                0029 #include "EEPARAMS.h"
                0030 #include "PARAMS.h"
                0031 #ifdef ALLOW_THSICE
                0032 #include "THSICE_PARAMS.h"
                0033 #include "THSICE_SIZE.h"
                0034 #include "THSICE_VARS.h"
                0035 #endif
                0036 #ifdef ALLOW_SEAICE
                0037 # include "SEAICE_SIZE.h"
                0038 # include "SEAICE.h"
                0039 #endif
                0040 
                0041       INTEGER siLo, siHi, sjLo, sjHi
                0042       PARAMETER ( siLo = 1-OLx , siHi = sNx+OLx )
                0043       PARAMETER ( sjLo = 1-OLy , sjHi = sNy+OLy )
                0044 
                0045 C     !INPUT PARAMETERS:
                0046 C     == Routine Arguments ==
                0047 C     swDown   :: incoming short-wave radiation (+=dw) [W/m2]
                0048 C     lwDown   :: incoming  long-wave radiation (+=dw) [W/m2]
                0049 C     uRelWind :: relative wind speed, u-component [m/s]
                0050 C     vRelWind :: relative wind speed, v-component [m/s]
                0051 C     LVapor   :: latent heat of vaporisation
                0052 C     bi, bj   :: tile indices
                0053 C     myIter   :: current iteration number
                0054 C     myTime   :: current time in simulation
                0055 C     myThid   :: my Thread Id number
                0056       _RL  swDown(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0057       _RL  lwDown(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0058       _RL uWind(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0059       _RL vWind(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0060       _RL LVapor
                0061       _RL myTime
                0062       INTEGER bi, bj, myIter, myThid
                0063 
                0064 C     !OUTPUT PARAMETERS:
                0065 C     fsha     :: sensible heat-flux over seaice (+=up) [W/m2]
                0066 C     flha     :: latent heat-flux over seaice   (+=up) [W/m2]
                0067 C     evp      :: evaporation over seaice     (+=up) [kg/m2/s]
                0068 C     xolw     :: upward long-wave over seaice   (+=up) [W/m2]
                0069 C     ssqt     ::
                0070 C     q100     ::
                0071 C     cdq      ::
                0072 C     Tsurf    :: updated seaice/snow surface temperature [deg.C]
                0073 C     iceFrac  :: ice fraction [0-1]
                0074 C     sw2oce   :: short-wave over seaice into the ocean (+=dw) [W/m2]
                0075       _RL fsha(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0076       _RL flha(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0077       _RL evp (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0078       _RL xolw(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0079       _RL ssqt(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0080       _RL q100(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0081       _RL cdq (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0082       _RL Tsurf(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0083       _RL iceFrac(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0084       _RS sw2oce (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0085 c     _RL prcAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0086 
                0087 #ifdef ALLOW_THSICE
                0088 C     !LOCAL VARIABLES:
                0089 C     == Local variables ==
                0090 C     uRelWind :: relative wind speed, u-component [m/s], (C-grid)
                0091 C     vRelWind :: relative wind speed, v-component [m/s], (C-grid)
                0092 C     windSq   :: relative wind speed squared (grid-cell center)
                0093       INTEGER i, j
                0094       INTEGER iceOrNot
                0095       INTEGER iMin, iMax
                0096       INTEGER jMin, jMax
                0097       _RL LatentHeat
                0098       _RL icFrac, opFrac
                0099         _RL netSW (1:sNx,1:sNy)
                0100         _RL sFlx  (1:sNx,1:sNy,0:2)
                0101 c       _RL tFrzOce(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0102         _RL dTsurf(1:sNx,1:sNy)
                0103 
                0104         _RL uRelWind(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0105         _RL vRelWind(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0106         _RL windSq(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0107         _RL cdu, dumArg(4)
                0108         _RL fsha0(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0109         _RL evp_0(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0110         _RL xolw0(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0111 c       _RL ssqt0(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0112 c       _RL q10_0(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0113 c       _RL cdq_0(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0114         _RL dShdTs(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0115         _RL dEvdTs(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0116         _RL dLwdTs(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0117 CEOP
                0118 
                0119 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0120 
                0121       iMin = 1
                0122       iMax = sNx
                0123       jMin = 1
                0124       jMax = sNy
                0125       LatentHeat = Lfresh + LVapor
                0126 
                0127 c     DO bj=myByLo(myThid),myByHi(myThid)
                0128 c      DO bi=myBxLo(myThid),myBxHi(myThid)
                0129 
                0130         CALL THSICE_GET_OCEAN(
                0131      I                         bi, bj, myTime, myIter, myThid )
                0132 
                0133         DO j = 1-OLy, sNy+OLy
                0134           DO i = 1-OLx, sNx+OLx
                0135             uRelWind(i,j) = uWind(i,j)
                0136             vRelWind(i,j) = vWind(i,j)
                0137           ENDDO
                0138         ENDDO
                0139 #ifdef ALLOW_SEAICE
                0140         IF ( useSEAICE ) THEN
                0141          DO j = 1-OLy, sNy+OLy
                0142           DO i = 1-OLx, sNx+OLx
                0143             uRelWind(i,j) = uRelWind(i,j)-uIce(i,j,bi,bj)
                0144             vRelWind(i,j) = vRelWind(i,j)-vIce(i,j,bi,bj)
                0145           ENDDO
                0146          ENDDO
                0147         ENDIF
                0148 #endif /* ALLOW_SEAICE */
                0149         DO j = jMin,jMax
                0150           DO i = iMin,iMax
                0151             windSq(i,j) = ( uRelWind( i ,j)*uRelWind( i ,j)
                0152      &                    + uRelWind(i+1,j)*uRelWind(i+1,j)
                0153      &                    + vRelWind(i, j )*vRelWind(i, j )
                0154      &                    + vRelWind(i,j+1)*vRelWind(i,j+1)
                0155      &                    )*0.5 _d 0
                0156           ENDDO
                0157         ENDDO
                0158 
                0159 C   1) compute albedo ; compute netSW
                0160        CALL THSICE_ALBEDO(
                0161      I          bi, bj, siLo, siHi, sjLo, sjHi,
                0162      I          iMin,iMax, jMin,jMax,
                0163      I          iceMask(siLo,sjLo,bi,bj), iceHeight(siLo,sjLo,bi,bj),
                0164      I          snowHeight(siLo,sjLo,bi,bj), Tsrf(siLo,sjLo,bi,bj),
                0165      I          snowAge(siLo,sjLo,bi,bj),
                0166      O          siceAlb(siLo,sjLo,bi,bj), icAlbNIR(siLo,sjLo,bi,bj),
                0167      I          myTime, myIter, myThid )
                0168 
                0169        DO j = jMin, jMax
                0170         DO i = iMin, iMax
                0171          IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN
                0172 C-      surface net SW flux:
                0173           netSW(i,j) = swDown(i,j)
                0174      &               *(1. _d 0 - siceAlb(i,j,bi,bj))
                0175          ELSE
                0176           netSW(i,j) = swDown(i,j)
                0177          ENDIF
                0178         ENDDO
                0179        ENDDO
                0180 
                0181 
                0182 C   2) compute other flx over seaice, over melting surf
                0183 C   3) compute other flx over seaice & derivative vs Tsurf, using previous Tsurf
                0184          DO j = jMin, jMax
                0185           DO i = iMin, iMax
                0186 
                0187             IF ( snowHeight(i,j,bi,bj).GT.3. _d -1 ) THEN
                0188              iceornot=2
                0189             ELSE
                0190              iceornot=1
                0191             ENDIF
                0192             Tsurf(i,j) = 0.
                0193             CALL CHEAPAML_COARE3_FLUX(
                0194      I                    i, j, bi, bj, iceOrNot,
                0195      I                    Tsurf, windSq,
                0196      O                    fsha0(i,j), flha(i,j), evp_0(i,j), xolw0(i,j),
                0197      O                    ssqt(i,j), q100(i,j), cdq(i,j), cdu,
                0198      O                    dumArg(1), dumArg(2), dumArg(3), dumArg(4),
                0199      I                    myIter, myThid )
                0200             sFlx(i,j,0) = lwDown(i,j)- xolw0(i,j)
                0201      &                  - fsha0(i,j) - evp_0(i,j)*LatentHeat
                0202 
                0203             Tsurf(i,j) = Tsrf(i,j,bi,bj)
                0204             CALL CHEAPAML_COARE3_FLUX(
                0205      I                    i, j, bi, bj, iceOrNot,
                0206      I                    Tsurf, windSq,
                0207      O                    fsha(i,j), flha(i,j), evp(i,j), xolw(i,j),
                0208      O                    ssqt(i,j), q100(i,j), cdq(i,j), cdu,
                0209      O              dShdTs(i,j), dEvdTs(i,j), dLwdTs(i,j), dumArg(4),
                0210      I                   myIter, myThid )
                0211             sFlx(i,j,1) = lwDown(i,j)- xolw(i,j)
                0212      &                  - fsha(i,j) - evp(i,j)*LatentHeat
                0213             sFlx(i,j,2) = -dLwdTs(i,j)
                0214      &                  - dShdTs(i,j) - dEvdTs(i,j)*LatentHeat
                0215           ENDDO
                0216          ENDDO
                0217 
                0218 C   4) solve for surf & seaice temp
                0219 C--    needs to fill in snowPrc, ( & prcAtm ? )
                0220 C--    note: this S/R  assumes No overlap
                0221          CALL THSICE_IMPL_TEMP(
                0222      I                netSW, sFlx,
                0223      O                dTsurf,
                0224      I                bi, bj, myTime, myIter, myThid )
                0225 
                0226 C   5) update surf fluxes
                0227         DO j = jMin, jMax
                0228          DO i = iMin, iMax
                0229           iceFrac(i,j) = iceMask(i,j,bi,bj)
                0230           sw2oce (i,j) = icFlxSW(i,j,bi,bj)
                0231           IF ( dTsurf(i,j) .GT. 999. ) THEN
                0232 c          dTsurf(J)= tFreeze - Tsurf(J)
                0233            Tsurf(i,j)= 0.
                0234            fsha(i,j) = fsha0(i,j)
                0235            flha(i,j) = evp_0(i,j)*LatentHeat
                0236            evp(i,j)  = evp_0(i,j)
                0237            xolw(i,j) = xolw0(i,j)
                0238           ELSE
                0239            Tsurf(i,j)= Tsurf(i,j)+ dTsurf(i,j)
                0240            fsha(i,j) = fsha(i,j) + dTsurf(i,j)*dShdTs(i,j)
                0241            evp(i,j)  = evp(i,j)  + dTsurf(i,j)*dEvdTs(i,j)
                0242            flha(i,j) = evp(i,j)*LatentHeat
                0243            xolw(i,j) = xolw(i,j) + dTsurf(i,j)*dLwdTs(i,j)
                0244           ENDIF
                0245          ENDDO
                0246         ENDDO
                0247 
                0248         DO j = jMin, jMax
                0249          DO i = iMin, iMax
                0250 c         IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN
                0251            icFrac  = iceMask(i,j,bi,bj)
                0252            opFrac = 1. _d 0 - icFrac
                0253 C--    Update Fluxes :
                0254            icFlxAtm(i,j,bi,bj) = netSW(i,j)
                0255      &                         + lwDown(i,j)- xolw(i,j)
                0256      &                         - fsha(i,j) - evp(i,j)*LVapor
                0257            icFrwAtm(i,j,bi,bj) = evp(i,j)
                0258 c         ENDIF
                0259          ENDDO
                0260         ENDDO
                0261 
                0262 c      ENDDO
                0263 c     ENDDO
                0264 
                0265 #endif /* ALLOW_THSICE */
                0266       RETURN
                0267       END