| 1 | dimitri | 1.1 | C $Header: /u/gcmpack/MITgcm/pkg/seaice/budget.F,v 1.11 2003/10/09 04:19:20 edhill Exp $ | 
| 2 |  |  | C $Name:  $ | 
| 3 |  |  |  | 
| 4 |  |  | #include "SEAICE_OPTIONS.h" | 
| 5 |  |  |  | 
| 6 |  |  | CStartOfInterface | 
| 7 |  |  | SUBROUTINE BUDGET(UG, TICE, HICE1, FICE1, KOPEN, bi, bj) | 
| 8 |  |  | C     /==========================================================\ | 
| 9 |  |  | C     | SUBROUTINE budget                                        | | 
| 10 |  |  | C     | o Calculate ice growth rate                              | | 
| 11 |  |  | C     |   see Hibler, MWR, 108, 1943-1973, 1980                  | | 
| 12 |  |  | C     |==========================================================| | 
| 13 |  |  | C     \==========================================================/ | 
| 14 |  |  | IMPLICIT NONE | 
| 15 |  |  |  | 
| 16 |  |  | C     === Global variables === | 
| 17 |  |  | #include "SIZE.h" | 
| 18 |  |  | #include "EEPARAMS.h" | 
| 19 |  |  | #include "FFIELDS.h" | 
| 20 |  |  | #include "SEAICE_PARAMS.h" | 
| 21 |  |  | #include "SEAICE_FFIELDS.h" | 
| 22 |  |  |  | 
| 23 |  |  | C     Subset of variables from SEAICE.h | 
| 24 |  |  | _RL HEFF       (1-OLx:sNx+OLx,1-OLy:sNy+OLy,3,nSx,nSy) | 
| 25 |  |  | _RL HSNOW      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,  nSx,nSy) | 
| 26 |  |  | _RL QNETO      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,  nSx,nSy) | 
| 27 |  |  | _RL QNETI      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,  nSx,nSy) | 
| 28 |  |  | _RL QSWO       (1-OLx:sNx+OLx,1-OLy:sNy+OLy,  nSx,nSy) | 
| 29 |  |  | _RL QSWI       (1-OLx:sNx+OLx,1-OLy:sNy+OLy,  nSx,nSy) | 
| 30 |  |  | COMMON/TRANS/HEFF,HSNOW | 
| 31 |  |  | COMMON/QFLUX/QNETO,QNETI,QSWO,QSWI | 
| 32 |  |  |  | 
| 33 |  |  | C     === Routine arguments === | 
| 34 |  |  | _RL UG         (1-OLx:sNx+OLx, 1-OLy:sNy+OLy) | 
| 35 |  |  | _RL TICE       (1-OLx:sNx+OLx, 1-OLy:sNy+OLy,  nSx,nSy) | 
| 36 |  |  | _RL HICE1      (1-OLx:sNx+OLx, 1-OLy:sNy+OLy) | 
| 37 |  |  | _RL FICE1       (1-OLx:sNx+OLx, 1-OLy:sNy+OLy,  nSx,nSy) | 
| 38 |  |  | INTEGER KOPEN | 
| 39 |  |  | INTEGER bi, bj | 
| 40 |  |  | CEndOfInterface | 
| 41 |  |  |  | 
| 42 |  |  | #ifdef ALLOW_SEAICE | 
| 43 |  |  |  | 
| 44 |  |  | C     === Local variables === | 
| 45 |  |  | C     i,j,k,bi,bj - Loop counters | 
| 46 |  |  |  | 
| 47 |  |  | INTEGER i, j | 
| 48 |  |  | INTEGER ITER | 
| 49 |  |  | _RL  QS1, C1, C2, C3, C4, C5, TB, D1, D1W, D1I, D3 | 
| 50 |  |  | _RL  TMELT, TMELTP, XKI, XKS, HCUT, ASNOW, XIO | 
| 51 |  |  |  | 
| 52 |  |  | _RL HICE       (1-OLx:sNx+OLx, 1-OLy:sNy+OLy) | 
| 53 |  |  | _RL ALB        (1-OLx:sNx+OLx, 1-OLy:sNy+OLy) | 
| 54 |  |  | _RL A1         (1-OLx:sNx+OLx, 1-OLy:sNy+OLy) | 
| 55 |  |  | _RL A2         (1-OLx:sNx+OLx, 1-OLy:sNy+OLy) | 
| 56 |  |  | _RL A3         (1-OLx:sNx+OLx, 1-OLy:sNy+OLy) | 
| 57 |  |  | _RL B          (1-OLx:sNx+OLx, 1-OLy:sNy+OLy) | 
| 58 |  |  |  | 
| 59 |  |  | _RL mymin_R8, mymax_R8 | 
| 60 |  |  | external mymin_R8, mymax_R8 | 
| 61 |  |  |  | 
| 62 |  |  | C IF KOPEN LT 0, THEN DO OPEN WATER BUDGET | 
| 63 |  |  | C NOW DEFINE ASSORTED CONSTANTS | 
| 64 |  |  | C SATURATION VAPOR PRESSURE CONSTANT | 
| 65 |  |  | QS1=0.622 _d +00/1013.0 _d +00 | 
| 66 |  |  | C MAYKUTS CONSTANTS FOR SAT. VAP. PRESSURE TEMP. POLYNOMIAL | 
| 67 |  |  | C1=2.7798202 _d -06 | 
| 68 |  |  | C2=-2.6913393 _d -03 | 
| 69 |  |  | C3=0.97920849 _d +00 | 
| 70 |  |  | C4=-158.63779 _d +00 | 
| 71 |  |  | C5=9653.1925 _d +00 | 
| 72 |  |  | C FREEZING TEMPERATURE OF SEAWATER | 
| 73 |  |  | TB=271.2 _d +00 | 
| 74 |  |  | C SENSIBLE HEAT CONSTANT | 
| 75 |  |  | D1=SEAICE_sensHeat | 
| 76 |  |  | C WATER LATENT HEAT CONSTANT | 
| 77 |  |  | D1W=SEAICE_latentWater | 
| 78 |  |  | C ICE LATENT HEAT CONSTANT | 
| 79 |  |  | D1I=SEAICE_latentIce | 
| 80 |  |  | C STEFAN BOLTZMAN CONSTANT TIMES 0.97 EMISSIVITY | 
| 81 |  |  | D3=SEAICE_emissivity | 
| 82 |  |  | C MELTING TEMPERATURE OF ICE | 
| 83 |  |  | TMELT=273.16 _d +00 | 
| 84 |  |  | TMELTP=273.159 _d +00 | 
| 85 |  |  | C ICE CONDUCTIVITY | 
| 86 |  |  | XKI=SEAICE_iceConduct | 
| 87 |  |  | C SNOW CONDUCTIVITY | 
| 88 |  |  | XKS=SEAICE_snowConduct | 
| 89 |  |  | C CUTOFF SNOW THICKNESS | 
| 90 |  |  | HCUT=SEAICE_snowThick | 
| 91 |  |  | C PENETRATION SHORTWAVE RADIATION FACTOR | 
| 92 |  |  | XIO=SEAICE_shortwave | 
| 93 |  |  |  | 
| 94 |  |  | DO J=1,sNy | 
| 95 |  |  | DO I=1,sNx | 
| 96 |  |  | TICE(I,J,bi,bj)=MYMIN_R8(273.16 _d 0+MAX_TICE,TICE(I,J,bi,bj)) | 
| 97 |  |  | ATEMP(I,J,bi,bj)=MYMAX_R8(273.16 _d 0+MIN_ATEMP,ATEMP(I,J,bi,bj)) | 
| 98 |  |  | LWDOWN(I,J,bi,bj)=MYMAX_R8(MIN_LWDOWN,LWDOWN(I,J,bi,bj)) | 
| 99 |  |  | ENDDO | 
| 100 |  |  | ENDDO | 
| 101 |  |  |  | 
| 102 |  |  | C NOW DECIDE IF OPEN WATER OR ICE | 
| 103 |  |  | IF(KOPEN.LE.0) THEN | 
| 104 |  |  |  | 
| 105 |  |  | C NOW DETERMINE OPEN WATER HEAT BUD. ASSUMING TICE=WATER TEMP. | 
| 106 |  |  | C WATER ALBEDO IS ASSUMED TO BE THE CONSTANT SEAICE_waterAlbedo | 
| 107 |  |  | DO J=1,sNy | 
| 108 |  |  | DO I=1,sNx | 
| 109 |  |  | #ifdef SEAICE_EXTERNAL_FLUXES | 
| 110 |  |  | FICE1(I,J,bi,bj)=QNET(I,J,bi,bj)+Qsw(I,J,bi,bj) | 
| 111 |  |  | QSWO(I,J,bi,bj)=Qsw(I,J,bi,bj) | 
| 112 |  |  | #else /* SEAICE_EXTERNAL_FLUXES undefined */ | 
| 113 |  |  | ALB(I,J)=SEAICE_waterAlbedo | 
| 114 |  |  | A1(I,J)=(ONE-ALB(I,J))*SWDOWN(I,J,bi,bj) | 
| 115 |  |  | &          +LWDOWN(I,J,bi,bj)*0.97 _d 0 | 
| 116 |  |  | &          +D1*UG(I,J)*ATEMP(I,J,bi,bj)+D1W*UG(I,J)*AQH(I,J,bi,bj) | 
| 117 |  |  | B(I,J)=QS1*6.11 _d +00*EXP(17.2694 _d +00 | 
| 118 |  |  | &          *(TICE(I,J,bi,bj)-TMELT) | 
| 119 |  |  | &          /(TICE(I,J,bi,bj)-TMELT+237.3 _d +00)) | 
| 120 |  |  | A2(I,J)=-D1*UG(I,J)*TICE(I,J,bi,bj)-D1W*UG(I,J)*B(I,J) | 
| 121 |  |  | &          -D3*(TICE(I,J,bi,bj)**4) | 
| 122 |  |  | FICE1(I,J,bi,bj)=-A1(I,J)-A2(I,J) | 
| 123 |  |  | QSWO(I,J,bi,bj)=-(ONE-ALB(I,J))*SWDOWN(I,J,bi,bj) | 
| 124 |  |  | #endif /* SEAICE_EXTERNAL_FLUXES */ | 
| 125 |  |  | QNETO(I,J,bi,bj)=FICE1(I,J,bi,bj)-QSWO(I,J,bi,bj) | 
| 126 |  |  | ENDDO | 
| 127 |  |  | ENDDO | 
| 128 |  |  |  | 
| 129 |  |  | ELSE | 
| 130 |  |  |  | 
| 131 |  |  | C COME HERE IF ICE COVER | 
| 132 |  |  | C FIRST PUT MINIMUM ON ICE THICKNESS | 
| 133 |  |  | DO J=1,sNy | 
| 134 |  |  | DO I=1,sNx | 
| 135 |  |  | HICE(I,J)=MYMAX_R8(HICE1(I,J),0.05 _d +00) | 
| 136 |  |  | HICE(I,J)=MYMIN_R8(HICE(I,J),9.0 _d +00) | 
| 137 |  |  | ENDDO | 
| 138 |  |  | ENDDO | 
| 139 |  |  | C NOW DECIDE ON ALBEDO | 
| 140 |  |  | DO J=1,sNy | 
| 141 |  |  | DO I=1,sNx | 
| 142 |  |  | ALB(I,J)=SEAICE_dryIceAlb | 
| 143 |  |  | IF(TICE(I,J,bi,bj).GT.TMELTP) ALB(I,J)=SEAICE_wetIceAlb | 
| 144 |  |  | ASNOW=SEAICE_drySnowAlb | 
| 145 |  |  | IF(TICE(I,J,bi,bj).GT.TMELTP) ASNOW=SEAICE_wetSnowAlb | 
| 146 |  |  | IF(HSNOW(I,J,bi,bj).GT.HCUT) THEN | 
| 147 |  |  | ALB(I,J)=ASNOW | 
| 148 |  |  | ELSE | 
| 149 |  |  | ALB(I,J)=ALB(I,J)+(HSNOW(I,J,bi,bj)/HCUT)*(ASNOW-ALB(I,J)) | 
| 150 |  |  | IF(ALB(I,J).GT.ASNOW) ALB(I,J)=ASNOW | 
| 151 |  |  | END IF | 
| 152 |  |  | ENDDO | 
| 153 |  |  | ENDDO | 
| 154 |  |  | C NOW DETERMINE FIXED FORCING TERM IN HEAT BUDGET | 
| 155 |  |  | DO J=1,sNy | 
| 156 |  |  | DO I=1,sNx | 
| 157 |  |  | IF(HSNOW(I,J,bi,bj).GT.0.0) THEN | 
| 158 |  |  | C NO SW PENETRATION WITH SNOW | 
| 159 |  |  | A1(I,J)=(ONE-ALB(I,J))*SWDOWN(I,J,bi,bj) | 
| 160 |  |  | &       +LWDOWN(I,J,bi,bj)*0.97 _d 0 | 
| 161 |  |  | &       +D1*UG(I,J)*ATEMP(I,J,bi,bj)+D1I*UG(I,J)*AQH(I,J,bi,bj) | 
| 162 |  |  | ELSE | 
| 163 |  |  | C SW PENETRATION UNDER ICE | 
| 164 |  |  | A1(I,J)=(ONE-ALB(I,J))*SWDOWN(I,J,bi,bj) | 
| 165 |  |  | &       *(ONE-XIO*EXP(-1.5 _d 0*HICE(I,J))) | 
| 166 |  |  | &       +LWDOWN(I,J,bi,bj)*0.97 _d 0 | 
| 167 |  |  | &       +D1*UG(I,J)*ATEMP(I,J,bi,bj)+D1I*UG(I,J)*AQH(I,J,bi,bj) | 
| 168 |  |  | ENDIF | 
| 169 |  |  | ENDDO | 
| 170 |  |  | ENDDO | 
| 171 |  |  | C NOW COMPUTE OTHER TERMS IN HEAT BUDGET | 
| 172 |  |  | C COME HERE AT START OF ITERATION | 
| 173 |  |  |  | 
| 174 |  |  | crg check wether a2 is needed in the list of variables | 
| 175 |  |  | cdm Ralf, the line below causes following error message | 
| 176 |  |  | cdm INTERNAL ERROR: cannot find var clone to ada2 | 
| 177 |  |  | cdm c$taf loop = iteration TICE,A2 | 
| 178 |  |  | cdm iterative solver for ice growth rate | 
| 179 |  |  | cdm inputs:  TICE  ice temperature | 
| 180 |  |  | cdm          UG    forcing | 
| 181 |  |  | cdm          HSNOW snow thickness | 
| 182 |  |  | cdm          HICE  ice thickness | 
| 183 |  |  | cdm outputs: A2 is needed for FICE1, which is ice growth rate | 
| 184 |  |  | cdm          TICE | 
| 185 |  |  | DO ITER=1,IMAX_TICE | 
| 186 |  |  |  | 
| 187 |  |  | DO J=1,sNy | 
| 188 |  |  | DO I=1,sNx | 
| 189 |  |  | B(I,J)=QS1*(C1*TICE(I,J,bi,bj)**4+C2*TICE(I,J,bi,bj)**3 | 
| 190 |  |  | &            +C3*TICE(I,J,bi,bj)**2+C4*TICE(I,J,bi,bj)+C5) | 
| 191 |  |  | A2(I,J)=-D1*UG(I,J)*TICE(I,J,bi,bj)-D1I*UG(I,J)*B(I,J) | 
| 192 |  |  | &             -D3*(TICE(I,J,bi,bj)**4) | 
| 193 |  |  | B(I,J)=XKS/(HSNOW(I,J,bi,bj)/HICE(I,J)+XKS/XKI)/HICE(I,J) | 
| 194 |  |  | A3(I,J)=4.0 _d +00*D3*(TICE(I,J,bi,bj)**3)+B(I,J)+D1*UG(I,J) | 
| 195 |  |  | B(I,J)=B(I,J)*(TB-TICE(I,J,bi,bj)) | 
| 196 |  |  | cdm | 
| 197 |  |  | cdm   if(TICE(I,J,bi,bj).le.206.) | 
| 198 |  |  | cdm  &           print '(A,3i4,f12.2)','### ITER,I,J,TICE', | 
| 199 |  |  | cdm  &           ITER,I,J,TICE(I,J,bi,bj) | 
| 200 |  |  | cdm | 
| 201 |  |  | ENDDO | 
| 202 |  |  | ENDDO | 
| 203 |  |  | C NOW DECIDE IF IT IS TIME TO ESTIMATE GROWTH RATES | 
| 204 |  |  | C NOW DETERMINE NEW ICE TEMPERATURE | 
| 205 |  |  | DO J=1,sNy | 
| 206 |  |  | DO I=1,sNx | 
| 207 |  |  | TICE(I,J,bi,bj)=TICE(I,J,bi,bj) | 
| 208 |  |  | &                     +(A1(I,J)+A2(I,J)+B(I,J))/A3(I,J) | 
| 209 |  |  | TICE(I,J,bi,bj)=MYMAX_R8(273.16 _d 0+MIN_TICE,TICE(I,J,bi,bj)) | 
| 210 |  |  | ENDDO | 
| 211 |  |  | ENDDO | 
| 212 |  |  | C NOW SET ICE TEMP TO MIN OF TMELT/ITERATION RESULT | 
| 213 |  |  | DO J=1,sNy | 
| 214 |  |  | DO I=1,sNx | 
| 215 |  |  | TICE(I,J,bi,bj)=MYMIN_R8(TICE(I,J,bi,bj),TMELT) | 
| 216 |  |  | ENDDO | 
| 217 |  |  | ENDDO | 
| 218 |  |  |  | 
| 219 |  |  | C END OF ITERATION | 
| 220 |  |  | ENDDO | 
| 221 |  |  |  | 
| 222 |  |  | DO J=1,sNy | 
| 223 |  |  | DO I=1,sNx | 
| 224 |  |  | FICE1(I,J,bi,bj)=-A1(I,J)-A2(I,J) | 
| 225 |  |  | IF(HSNOW(I,J,bi,bj).GT.0.0) THEN | 
| 226 |  |  | C NO SW PENETRATION WITH SNOW | 
| 227 |  |  | QSWI(I,J,bi,bj)=ZERO | 
| 228 |  |  | ELSE | 
| 229 |  |  | C SW PENETRATION UNDER ICE | 
| 230 |  |  | QSWI(I,J,bi,bj)=-(ONE-ALB(I,J))*SWDOWN(I,J,bi,bj) | 
| 231 |  |  | &       *XIO*EXP(-1.5 _d 0*HICE(I,J)) | 
| 232 |  |  | ENDIF | 
| 233 |  |  | ENDDO | 
| 234 |  |  | ENDDO | 
| 235 |  |  |  | 
| 236 |  |  | END IF | 
| 237 |  |  |  | 
| 238 |  |  | #endif /* ALLOW_SEAICE */ | 
| 239 |  |  |  | 
| 240 |  |  | RETURN | 
| 241 |  |  | END |