| 1 | C $Header: /u/gcmpack/MITgcm/pkg/seaice/growth.F,v 1.19 2004/07/06 06:35:04 dimitri Exp $ | 
| 2 | C $Name:  $ | 
| 3 |  | 
| 4 | #include "SEAICE_OPTIONS.h" | 
| 5 |  | 
| 6 | CStartOfInterface | 
| 7 | SUBROUTINE growth( myTime, myIter, myThid ) | 
| 8 | C     /==========================================================\ | 
| 9 | C     | SUBROUTINE growth                                        | | 
| 10 | C     | o Updata ice thickness and snow depth                    | | 
| 11 | C     |==========================================================| | 
| 12 | C     \==========================================================/ | 
| 13 | IMPLICIT NONE | 
| 14 |  | 
| 15 | C     === Global variables === | 
| 16 | #include "SIZE.h" | 
| 17 | #include "EEPARAMS.h" | 
| 18 | #include "PARAMS.h" | 
| 19 | #include "DYNVARS.h" | 
| 20 | #include "GRID.h" | 
| 21 | #include "FFIELDS.h" | 
| 22 | #include "SEAICE_PARAMS.h" | 
| 23 | #include "SEAICE.h" | 
| 24 | #include "SEAICE_FFIELDS.h" | 
| 25 |  | 
| 26 | #ifdef ALLOW_AUTODIFF_TAMC | 
| 27 | # include "tamc.h" | 
| 28 | #endif | 
| 29 | C     === Routine arguments === | 
| 30 | C     myTime - Simulation time | 
| 31 | C     myIter - Simulation timestep number | 
| 32 | C     myThid - Thread no. that called this routine. | 
| 33 | _RL myTime | 
| 34 | INTEGER myIter, myThid | 
| 35 | CEndOfInterface | 
| 36 |  | 
| 37 | #ifdef ALLOW_SEAICE | 
| 38 |  | 
| 39 | C     === Local variables === | 
| 40 | C     i,j,bi,bj - Loop counters | 
| 41 |  | 
| 42 | INTEGER i, j, bi, bj | 
| 43 | _RL  TBC, salinity_ice, SDF, Q0, QS | 
| 44 | _RL GAREA( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy           ) | 
| 45 | _RL GHEFF( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy           ) | 
| 46 | _RL AR   ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, nSx, nSy ) | 
| 47 |  | 
| 48 | C     number of surface interface layer | 
| 49 | INTEGER kSurface | 
| 50 |  | 
| 51 | if ( buoyancyRelation .eq. 'OCEANICP' ) then | 
| 52 | kSurface        = Nr | 
| 53 | else | 
| 54 | kSurface        = 1 | 
| 55 | endif | 
| 56 |  | 
| 57 | salinity_ice=4.0 _d 0      ! ICE SALINITY | 
| 58 | TBC=271.2 _d 0-273.16 _d 0 ! FREEZING TEMP. OF SEA WATER | 
| 59 | SDF=1000.0 _d 0/330.0 _d 0 ! RATIO OF WATER DESITY AND SNOW DENSITY | 
| 60 | Q0=1.0D-06/302.0 _d +00    ! INVERSE HEAT OF FUSION OF ICE | 
| 61 | QS=1.1 _d +08              ! HEAT OF FUSION OF SNOW | 
| 62 |  | 
| 63 | DO bj=myByLo(myThid),myByHi(myThid) | 
| 64 | DO bi=myBxLo(myThid),myBxHi(myThid) | 
| 65 | c | 
| 66 | cph( | 
| 67 | #ifdef ALLOW_AUTODIFF_TAMC | 
| 68 | act1 = bi - myBxLo(myThid) | 
| 69 | max1 = myBxHi(myThid) - myBxLo(myThid) + 1 | 
| 70 | act2 = bj - myByLo(myThid) | 
| 71 | max2 = myByHi(myThid) - myByLo(myThid) + 1 | 
| 72 | act3 = myThid - 1 | 
| 73 | max3 = nTx*nTy | 
| 74 | act4 = ikey_dynamics - 1 | 
| 75 | iicekey = (act1 + 1) + act2*max1 | 
| 76 | &                      + act3*max1*max2 | 
| 77 | &                      + act4*max1*max2*max3 | 
| 78 | #endif /* ALLOW_AUTODIFF_TAMC */ | 
| 79 | c | 
| 80 | #ifdef ALLOW_AUTODIFF_TAMC | 
| 81 | CADJ STORE theta(:,:,:,bi,bj)= comlev1_bibj, | 
| 82 | CADJ &                         key = iicekey, byte = isbyte | 
| 83 | CADJ STORE area(:,:,:,bi,bj) = comlev1_bibj, | 
| 84 | CADJ &                         key = iicekey, byte = isbyte | 
| 85 | CADJ STORE atemp(:,:,bi,bj)  = comlev1_bibj, | 
| 86 | CADJ &                         key = iicekey, byte = isbyte | 
| 87 | #endif /* ALLOW_AUTODIFF_TAMC */ | 
| 88 | cph) | 
| 89 | DO J=1,sNy | 
| 90 | DO I=1,sNx | 
| 91 | SEAICE_SALT(I,J,bi,bj)=ZERO | 
| 92 | ENDDO | 
| 93 | ENDDO | 
| 94 | #ifdef ALLOW_AUTODIFF_TAMC | 
| 95 | CADJ STORE area(:,:,:,bi,bj) = comlev1_bibj, | 
| 96 | CADJ &                         key = iicekey, byte = isbyte | 
| 97 | CADJ STORE heff(:,:,:,bi,bj) = comlev1_bibj, | 
| 98 | CADJ &                         key = iicekey, byte = isbyte | 
| 99 | #endif /* ALLOW_AUTODIFF_TAMC */ | 
| 100 | DO J=1,sNy | 
| 101 | DO I=1,sNx | 
| 102 | AR(I,J,bi,bj)=MIN(AREA(I,J,2,bi,bj), | 
| 103 | &         HEFF(I,J,2,bi,bj)*1.0 _d +04) | 
| 104 | ENDDO | 
| 105 | ENDDO | 
| 106 | #ifdef ALLOW_AUTODIFF_TAMC | 
| 107 | CADJ STORE heff(:,:,:,bi,bj) = comlev1_bibj, | 
| 108 | CADJ &                         key = iicekey, byte = isbyte | 
| 109 | #endif /* ALLOW_AUTODIFF_TAMC */ | 
| 110 | DO J=1,sNy | 
| 111 | DO I=1,sNx | 
| 112 | C--   Create or melt sea-ice so that first-level oceanic temperature | 
| 113 | C     is approximately at the freezing point when there is sea-ice. | 
| 114 | C     Initially the units of YNEG are m of sea-ice. | 
| 115 | C     The factor dRf(1)/72.0764, used to convert temperature | 
| 116 | C     change in deg K to m of sea-ice, is approximately: | 
| 117 | C     dRf(1) * (sea water heat capacity = 3996 J/kg/K) | 
| 118 | C        * (density of sea-water = 1026 kg/m^3) | 
| 119 | C        / (latent heat of fusion of sea-ice = 334000 J/kg) | 
| 120 | C        / (density of sea-ice = 910 kg/m^3) | 
| 121 | C     Negative YNEG leads to ice growth. | 
| 122 | C     Positive YNEG leads to ice melting. | 
| 123 | YNEG(I,J,bi,bj)=(theta(I,J,1,bi,bj)-TBC)*.01 | 
| 124 | &           *dRf(1)/72.0764 _d 0 | 
| 125 | GHEFF(I,J)=HEFF(I,J,1,bi,bj) | 
| 126 | HEFF(I,J,1,bi,bj)=MAX(ZERO,HEFF(I,J,1,bi,bj)-YNEG(I,J,bi,bj)) | 
| 127 | YNEG(I,J,bi,bj)=GHEFF(I,J)-HEFF(I,J,1,bi,bj) | 
| 128 | SEAICE_SALT(I,J,bi,bj)=SEAICE_SALT(I,J,bi,bj)-YNEG(I,J,bi,bj) | 
| 129 | C--   Now convert YNEG back to deg K. | 
| 130 | YNEG(I,J,bi,bj)=YNEG(I,J,bi,bj)*recip_dRf(1)*72.0764 _d 0 | 
| 131 | ENDDO | 
| 132 | ENDDO | 
| 133 | c | 
| 134 | ENDDO | 
| 135 | ENDDO | 
| 136 |  | 
| 137 | cph( | 
| 138 | #ifdef ALLOW_AUTODIFF_TAMC | 
| 139 | CADJ STORE area   = comlev1, key = ikey_dynamics | 
| 140 | CADJ STORE atemp  = comlev1, key = ikey_dynamics | 
| 141 | CADJ STORE heff   = comlev1, key = ikey_dynamics | 
| 142 | CADJ STORE hsnow  = comlev1, key = ikey_dynamics | 
| 143 | CADJ STORE lwdown = comlev1, key = ikey_dynamics | 
| 144 | CADJ STORE tice   = comlev1, key = ikey_dynamics | 
| 145 | CADJ STORE uwind  = comlev1, key = ikey_dynamics | 
| 146 | CADJ STORE vwind  = comlev1, key = ikey_dynamics | 
| 147 | # ifdef SEAICE_MULTILEVEL | 
| 148 | CADJ STORE tices  = comlev1, key = ikey_dynamics | 
| 149 | # endif | 
| 150 | #endif /* ALLOW_AUTODIFF_TAMC */ | 
| 151 | cph) | 
| 152 | C GROWTH SUBROUTINE CALCULATES TOTAL GROWTH TENDENCIES, | 
| 153 | C INCLUDING SNOWFALL | 
| 154 | CALL GROATB(A22,myThid) | 
| 155 |  | 
| 156 | DO bj=myByLo(myThid),myByHi(myThid) | 
| 157 | DO bi=myBxLo(myThid),myBxHi(myThid) | 
| 158 | cph( | 
| 159 | #ifdef ALLOW_AUTODIFF_TAMC | 
| 160 | act1 = bi - myBxLo(myThid) | 
| 161 | max1 = myBxHi(myThid) - myBxLo(myThid) + 1 | 
| 162 | act2 = bj - myByLo(myThid) | 
| 163 | max2 = myByHi(myThid) - myByLo(myThid) + 1 | 
| 164 | act3 = myThid - 1 | 
| 165 | max3 = nTx*nTy | 
| 166 | act4 = ikey_dynamics - 1 | 
| 167 | iicekey = (act1 + 1) + act2*max1 | 
| 168 | &                      + act3*max1*max2 | 
| 169 | &                      + act4*max1*max2*max3 | 
| 170 | #endif /* ALLOW_AUTODIFF_TAMC */ | 
| 171 | c | 
| 172 | #ifdef ALLOW_AUTODIFF_TAMC | 
| 173 | CADJ STORE area(:,:,:,bi,bj) = comlev1_bibj, | 
| 174 | CADJ &                         key = iicekey, byte = isbyte | 
| 175 | CADJ STORE heff(:,:,:,bi,bj) = comlev1_bibj, | 
| 176 | CADJ &                         key = iicekey, byte = isbyte | 
| 177 | CADJ STORE hsnow(:,:,bi,bj)  = comlev1_bibj, | 
| 178 | CADJ &                         key = iicekey, byte = isbyte | 
| 179 | CADJ STORE fo(:,:,bi,bj)     = comlev1_bibj, | 
| 180 | CADJ &                         key = iicekey, byte = isbyte | 
| 181 | CADJ STORE fice(:,:,bi,bj)   = comlev1_bibj, | 
| 182 | CADJ &                         key = iicekey, byte = isbyte | 
| 183 | #endif /* ALLOW_AUTODIFF_TAMC */ | 
| 184 | cph) | 
| 185 | C NOW CALCULATE CORRECTED GROWTH | 
| 186 | DO J=1,sNy | 
| 187 | DO I=1,sNx | 
| 188 | GHEFF(I,J)=-SEAICE_deltaTtherm*FICE(I,J,bi,bj) | 
| 189 | GAREA(I,J)=HSNOW(I,J,bi,bj)*QS | 
| 190 | IF(GHEFF(I,J).GT.ZERO.AND.GHEFF(I,J).LE.GAREA(I,J)) THEN | 
| 191 | HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)-GHEFF(I,J)/QS | 
| 192 | C SNOW CONVERTED INTO WATER AND THEN INTO ICE | 
| 193 | C The factor 0.920 is used to convert m of sea-ice to m of freshwater | 
| 194 | SEAICE_SALT(I,J,bi,bj)=SEAICE_SALT(I,J,bi,bj) | 
| 195 | &                  -(GHEFF(I,J)/QS)/SDF/0.920 _d 0*AR(I,J,bi,bj) | 
| 196 | FICE(I,J,bi,bj)=ZERO | 
| 197 | ELSE IF(GHEFF(I,J).GT.GAREA(I,J)) THEN | 
| 198 | FICE(I,J,bi,bj)=-(GHEFF(I,J)-GAREA(I,J))/SEAICE_deltaTtherm | 
| 199 | SEAICE_SALT(I,J,bi,bj)=SEAICE_SALT(I,J,bi,bj) | 
| 200 | &               -HSNOW(I,J,bi,bj)/SDF/0.920 _d 0*AR(I,J,bi,bj) | 
| 201 | HSNOW(I,J,bi,bj)=0.0 | 
| 202 | END IF | 
| 203 |  | 
| 204 | ENDDO | 
| 205 | ENDDO | 
| 206 |  | 
| 207 | C NOW GET TOTAL GROWTH RATE | 
| 208 | DO J=1,sNy | 
| 209 | DO I=1,sNx | 
| 210 | FHEFF(I,J,bi,bj)=FICE(I,J,bi,bj)*AR(I,J,bi,bj) | 
| 211 | &                    +(ONE-AR(I,J,bi,bj))*FO(I,J,bi,bj) | 
| 212 | ENDDO | 
| 213 | ENDDO | 
| 214 |  | 
| 215 |  | 
| 216 | C NOW UPDATE AREA | 
| 217 | DO J=1,sNy | 
| 218 | DO I=1,sNx | 
| 219 | GHEFF(I,J)=-SEAICE_deltaTtherm*FHEFF(I,J,bi,bj)*Q0 | 
| 220 | GAREA(I,J)=SEAICE_deltaTtherm*FO(I,J,bi,bj)*Q0 | 
| 221 | GHEFF(I,J)=-ONE*MIN(HEFF(I,J,1,bi,bj),GHEFF(I,J)) | 
| 222 | GAREA(I,J)=MAX(ZERO,GAREA(I,J)) | 
| 223 | HCORR(I,J,bi,bj)=MIN(ZERO,GHEFF(I,J)) | 
| 224 | ENDDO | 
| 225 | ENDDO | 
| 226 | DO J=1,sNy | 
| 227 | DO I=1,sNx | 
| 228 | GAREA(I,J)=TWO*(ONE-AREA(I,J,2,bi,bj))*GAREA(I,J)/HO | 
| 229 | &    +HALF*HCORR(I,J,bi,bj)*AREA(I,J,2,bi,bj) | 
| 230 | &    /(HEFF(I,J,1,bi,bj)+.00001 _d 0) | 
| 231 | AREA(I,J,1,bi,bj)=AREA(I,J,1,bi,bj)+GAREA(I,J) | 
| 232 | ENDDO | 
| 233 | ENDDO | 
| 234 |  | 
| 235 | C NOW UPDATE HEFF | 
| 236 | DO J=1,sNy | 
| 237 | DO I=1,sNx | 
| 238 | GHEFF(I,J)=-SEAICE_deltaTtherm* | 
| 239 | &               FICE(I,J,bi,bj)*Q0*AR(I,J,bi,bj) | 
| 240 | GHEFF(I,J)=-ONE*MIN(HEFF(I,J,1,bi,bj),GHEFF(I,J)) | 
| 241 | HEFF(I,J,1,bi,bj)=HEFF(I,J,1,bi,bj)+GHEFF(I,J) | 
| 242 | SEAICE_SALT(I,J,bi,bj)=SEAICE_SALT(I,J,bi,bj)+GHEFF(I,J) | 
| 243 | C NOW CALCULATE QNETI UNDER ICE IF ANY | 
| 244 | QNETI(I,J,bi,bj)=(GHEFF(I,J)-SEAICE_deltaTtherm* | 
| 245 | &         FICE(I,J,bi,bj)*Q0*AR(I,J,bi,bj))/Q0/SEAICE_deltaTtherm | 
| 246 | ENDDO | 
| 247 | ENDDO | 
| 248 |  | 
| 249 | C NOW GET TOTAL QNET AND QSW | 
| 250 | DO J=1,sNy | 
| 251 | DO I=1,sNx | 
| 252 | QNET(I,J,bi,bj)=QNETI(I,J,bi,bj)*AR(I,J,bi,bj) | 
| 253 | &                    +(ONE-AR(I,J,bi,bj))*QNETO(I,J,bi,bj) | 
| 254 | QSW(I,J,bi,bj)=QSWI(I,J,bi,bj)*AR(I,J,bi,bj) | 
| 255 | &                    +(ONE-AR(I,J,bi,bj))*QSWO(I,J,bi,bj) | 
| 256 | #ifndef SHORTWAVE_HEATING | 
| 257 | QNET(I,J,bi,bj)=QNET(I,J,bi,bj)+QSW(I,J,bi,bj) | 
| 258 | #endif | 
| 259 | C Add YNEG contribution to QNET | 
| 260 | QNET(I,J,bi,bj)=QNET(I,J,bi,bj) | 
| 261 | &         +YNEG(I,J,bi,bj)/SEAICE_deltaTtherm*maskC(I,J,1,bi,bj) | 
| 262 | &         *HeatCapacity_Cp*recip_horiVertRatio*rhoConst | 
| 263 | &         *drF(kSurface)*hFacC(i,j,kSurface,bi,bj) | 
| 264 | ENDDO | 
| 265 | ENDDO | 
| 266 |  | 
| 267 | C NOW UPDATE OTHER THINGS | 
| 268 | DO J=1,sNy | 
| 269 | DO I=1,sNx | 
| 270 | IF(FICE(I,J,bi,bj).GT.ZERO) THEN | 
| 271 | C FREEZING, PRECIP ADDED AS SNOW | 
| 272 | HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)+SEAICE_deltaTtherm* | 
| 273 | &            PRECIP(I,J,bi,bj)*AREA(I,J,2,bi,bj)*SDF | 
| 274 | ELSE | 
| 275 | C ADD PRECIP AS RAIN, WATER CONVERTED INTO ICE BY /0.920 _d 0 | 
| 276 | SEAICE_SALT(I,J,bi,bj)=SEAICE_SALT(I,J,bi,bj) | 
| 277 | &            -PRECIP(I,J,bi,bj)*AREA(I,J,2,bi,bj)* | 
| 278 | &            SEAICE_deltaTtherm/0.920 _d 0 | 
| 279 | ENDIF | 
| 280 | c Now add in precip over open water directly into ocean as negative salt | 
| 281 | SEAICE_SALT(I,J,bi,bj)=SEAICE_SALT(I,J,bi,bj) | 
| 282 | &         -PRECIP(I,J,bi,bj)*(ONE-AREA(I,J,2,bi,bj)) | 
| 283 | &         *SEAICE_deltaTtherm/0.920 _d 0 | 
| 284 | C NOW GET FRESH WATER FLUX | 
| 285 | EmPmR(I,J,bi,bj)= maskC(I,J,1,bi,bj)*( | 
| 286 | &         EVAP(I,J,bi,bj)-RUNOFF(I,J,bi,bj) | 
| 287 | &         +SEAICE_SALT(I,J,bi,bj)*0.92 _d 0/SEAICE_deltaTtherm | 
| 288 | &         ) | 
| 289 | ENDDO | 
| 290 | ENDDO | 
| 291 |  | 
| 292 | #ifdef SEAICE_DEBUG | 
| 293 | c      CALL PLOT_FIELD_XYRS( UWIND,'Current UWIND ', myIter, myThid ) | 
| 294 | c      CALL PLOT_FIELD_XYRS( VWIND,'Current VWIND ', myIter, myThid ) | 
| 295 | CALL PLOT_FIELD_XYRS( GWATX,'Current GWATX ', myIter, myThid ) | 
| 296 | CALL PLOT_FIELD_XYRS( GWATY,'Current GWATY ', myIter, myThid ) | 
| 297 | CALL PLOT_FIELD_XYRL( FO,'Current FO ', myIter, myThid ) | 
| 298 | CALL PLOT_FIELD_XYRL( FHEFF,'Current FHEFF ', myIter, myThid ) | 
| 299 | CALL PLOT_FIELD_XYRL( QSW,'Current QSW ', myIter, myThid ) | 
| 300 | CALL PLOT_FIELD_XYRL( QNET,'Current QNET ', myIter, myThid ) | 
| 301 | CALL PLOT_FIELD_XYRL( EmPmR,'Current EmPmR ', myIter, myThid ) | 
| 302 | DO j=1-OLy,sNy+OLy | 
| 303 | DO i=1-OLx,sNx+OLx | 
| 304 | GHEFF(I,J)=SQRT(UICE(I,J,1,bi,bj)**2+VICE(I,J,1,bi,bj)**2) | 
| 305 | GAREA(I,J)=HEFF(I,J,1,bi,bj) | 
| 306 | print*,'I J QNET:',I, J, QNET(i,j,bi,bj), QSW(I,J,bi,bj) | 
| 307 | ENDDO | 
| 308 | ENDDO | 
| 309 | CALL PLOT_FIELD_XYRL( GHEFF,'Current UICE ', myIter, myThid ) | 
| 310 | CALL PLOT_FIELD_XYRL( GAREA,'Current HEFF ', myIter, myThid ) | 
| 311 | DO j=1-OLy,sNy+OLy | 
| 312 | DO i=1-OLx,sNx+OLx | 
| 313 | if(HEFF(i,j,1,bi,bj).gt.1.) then | 
| 314 | print '(A,2i4,3f10.2)','#### i j heff theta yneg',i,j, | 
| 315 | &            HEFF(i,j,1,bi,bj),theta(I,J,1,bi,bj),yneg(I,J,bi,bj) | 
| 316 | print '(A,3f10.2)','QSW, QNET before/after correction', | 
| 317 | &            QSW(I,J,bi,bj),QNETI(I,J,bi,bj)*AR(I,J,bi,bj) | 
| 318 | &           +(ONE-AR(I,J,bi,bj))*QNETO(I,J,bi,bj), QNET(I,J,bi,bj) | 
| 319 | endif | 
| 320 | ENDDO | 
| 321 | ENDDO | 
| 322 | #endif /* SEAICE_DEBUG */ | 
| 323 |  | 
| 324 | crg Added by Ralf Giering: do we need DO_WE_NEED_THIS ? | 
| 325 | #define DO_WE_NEED_THIS | 
| 326 | C NOW ZERO OUTSIDE POINTS | 
| 327 | DO J=1,sNy | 
| 328 | DO I=1,sNx | 
| 329 | C NOW SET AREA(I,J,1,bi,bj)=0 WHERE NO ICE IS | 
| 330 | AREA(I,J,1,bi,bj)=MIN(AREA(I,J,1,bi,bj) | 
| 331 | &                         ,HEFF(I,J,1,bi,bj)/.0001 _d 0) | 
| 332 | ENDDO | 
| 333 | ENDDO | 
| 334 | #ifdef ALLOW_AUTODIFF_TAMC | 
| 335 | CADJ STORE area(:,:,:,bi,bj) = comlev1_bibj, | 
| 336 | CADJ &                         key = iicekey, byte = isbyte | 
| 337 | #endif /* ALLOW_AUTODIFF_TAMC */ | 
| 338 | DO J=1,sNy | 
| 339 | DO I=1,sNx | 
| 340 | C NOW TRUNCATE AREA | 
| 341 | #ifdef DO_WE_NEED_THIS | 
| 342 | AREA(I,J,1,bi,bj)=MIN(ONE,AREA(I,J,1,bi,bj)) | 
| 343 | ENDDO | 
| 344 | ENDDO | 
| 345 | #ifdef ALLOW_AUTODIFF_TAMC | 
| 346 | CADJ STORE area(:,:,:,bi,bj) = comlev1_bibj, | 
| 347 | CADJ &                         key = iicekey, byte = isbyte | 
| 348 | #endif /* ALLOW_AUTODIFF_TAMC */ | 
| 349 | DO J=1,sNy | 
| 350 | DO I=1,sNx | 
| 351 | AREA(I,J,1,bi,bj)=MAX(ZERO,AREA(I,J,1,bi,bj)) | 
| 352 | HSNOW(I,J,bi,bj)=MAX(ZERO,HSNOW(I,J,bi,bj)) | 
| 353 | #endif | 
| 354 | AREA(I,J,1,bi,bj)=AREA(I,J,1,bi,bj)*HEFFM(I,J,bi,bj) | 
| 355 | HEFF(I,J,1,bi,bj)=HEFF(I,J,1,bi,bj)*HEFFM(I,J,bi,bj) | 
| 356 | #ifdef DO_WE_NEED_THIS | 
| 357 | c          HEFF(I,J,1,bi,bj)=MIN(MAX_HEFF,HEFF(I,J,1,bi,bj)) | 
| 358 | #endif | 
| 359 | HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)*HEFFM(I,J,bi,bj) | 
| 360 | ENDDO | 
| 361 | ENDDO | 
| 362 |  | 
| 363 | ENDDO | 
| 364 | ENDDO | 
| 365 |  | 
| 366 | #endif /* ALLOW_SEAICE */ | 
| 367 |  | 
| 368 | RETURN | 
| 369 | END |