| 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 |