| 1 | dimitri | 1.2 | C $Header: /u/gcmpack/MITgcm_contrib/lab_sea_test/growth.F,v 1.1 2004/07/11 06:19:16 dimitri Exp $ | 
| 2 | dimitri | 1.1 | 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 | dimitri | 1.2 | _RL mymin_R8, mymax_R8 | 
| 52 |  |  | external mymin_R8, mymax_R8 | 
| 53 |  |  |  | 
| 54 | dimitri | 1.1 | 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 | dimitri | 1.2 | AR(I,J,bi,bj)=MYMIN_R8(AREA(I,J,2,bi,bj), | 
| 106 | dimitri | 1.1 | &         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 | dimitri | 1.2 | HEFF(I,J,1,bi,bj)=MYMAX_R8(ZERO,HEFF(I,J,1,bi,bj)-YNEG(I,J,bi,bj)) | 
| 130 | dimitri | 1.1 | 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 | dimitri | 1.2 | 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 | dimitri | 1.1 | 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 | dimitri | 1.2 | GHEFF(I,J)=-ONE*MYMIN_R8(HEFF(I,J,1,bi,bj),GHEFF(I,J)) | 
| 244 | dimitri | 1.1 | 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 | dimitri | 1.2 | AREA(I,J,1,bi,bj)=MYMIN_R8(AREA(I,J,1,bi,bj) | 
| 334 | dimitri | 1.1 | &                         ,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 | dimitri | 1.2 | AREA(I,J,1,bi,bj)=MYMIN_R8(ONE,AREA(I,J,1,bi,bj)) | 
| 346 | dimitri | 1.1 | 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 | dimitri | 1.2 | 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 | dimitri | 1.1 | #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 | dimitri | 1.2 | c          HEFF(I,J,1,bi,bj)=MYMIN_R8(MAX_HEFF,HEFF(I,J,1,bi,bj)) | 
| 361 | dimitri | 1.1 | #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 |