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