File indexing completed on 2024-06-04 05:11:06 UTC
view on githubraw file Latest commit 54f5c8e3 on 2024-06-03 18:10:07 UTC
00f81e6785 Ou W*0001 #include "STIC_OPTIONS.h"
0002 #ifdef ALLOW_SHELFICE
0003 # include "SHELFICE_OPTIONS.h"
0004 #endif
0005 #ifdef ALLOW_AUTODIFF
0006 # include "AUTODIFF_OPTIONS.h"
0007 #endif
0008 #ifdef ALLOW_CTRL
0009 # include "CTRL_OPTIONS.h"
0010 #endif
0011
0012
0013
0014
0015 SUBROUTINE STIC_THERMODYNAMICS(
0016 I myTime, myIter, myThid )
0017
0018
0019
0020
0021
0022
0023
0024
0025
54f5c8e380 Jean*0026
00f81e6785 Ou W*0027
0028
0029
0030
0031
0032
0033 IMPLICIT NONE
0034
0035
0036 #include "SIZE.h"
0037 #include "EEPARAMS.h"
0038 #include "PARAMS.h"
0039 #include "EOS.h"
0040 #include "GRID.h"
0041 #include "DYNVARS.h"
0042 #include "FFIELDS.h"
0043 #ifdef ALLOW_SHELFICE
0044 # include "SHELFICE.h"
0045 # include "SHELFICE_COST.h"
0046 #endif
0047 #include "STIC.h"
0048 #ifdef ALLOW_CTRL
0049 # include "CTRL_SIZE.h"
0050 # include "CTRL.h"
0051 # include "CTRL_DUMMY.h"
0052 # ifdef ALLOW_GENTIM2D_CONTROL
0053 # include "CTRL_GENARR.h"
0054 # endif
0055 #endif /* ALLOW_CTRL */
0056 #ifdef ALLOW_AUTODIFF_TAMC
0057 # include "tamc.h"
0058 #endif /* ALLOW_AUTODIFF_TAMC */
0059
0060
0061
0062
54f5c8e380 Jean*0063
00f81e6785 Ou W*0064
0065 _RL myTime
0066 INTEGER myIter
0067 INTEGER myThid
0068
0069
0070
0071
0072
0073
0074
0075
0076
0077
0078
0079
0080
0081
0082
54f5c8e380 Jean*0083 INTEGER i, j, k
0084 INTEGER bi, bj
00f81e6785 Ou W*0085 INTEGER CURI, CURJ, FRONT_K
0086
0087 _RL tLoc
0088 _RL sLoc
0089 _RL pLoc
0090 _RL gammaTLoc
0091 _RL gammaSLoc
0092
0093 _RL ice_bottom_Z_C
0094 _RL wet_top_Z_N, wet_bottom_Z_N
0095 _RL iceFrontWetContact_Z_max
0096 _RL iceFrontVertContactFrac, iceFrontCellThickness
0097 _RL iceFrontWidth, iceFrontFaceArea
0098 _RL thermalConductionDistance, thermalConductionTemp
0099 _RL tmpHeatFlux, tmpFWFLX
0100 _RL tmpForcingT, tmpForcingS
0101 _RL tmpFac, icfgridareaFrac
0102 _RL tmpHeatFluxscaled, tmpFWFLXscaled
0103 INTEGER SI
0104 _RL insituT
0105
0106 #ifdef ALLOW_DIAGNOSTICS
0107 # ifdef SHI_ALLOW_GAMMAFRICT
0108 _RL uStarDiag(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
0109 # endif
0110 #ifdef NONLIN_FRSURF
0111 _RL tmpDiag(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
0112 _RL tmpDiagT(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
0113 #endif
0114 _RL tmpDiagIcfForcingT(1-OLx:sNx+OLx,
0115 & 1-OLy:sNy+OLy,Nr,nSx,nSy)
0116 _RL tmpDiagIcfForcingS(1-OLx:sNx+OLx,
0117 & 1-OLy:sNy+OLy,Nr,nSx,nSy)
0118 _RL tmpDiagShelficeForcingT(1-OLx:sNx+OLx,
0119 & 1-OLy:sNy+OLy,nSx,nSy)
0120 _RL tmpDiagShelficeForcingS(1-OLx:sNx+OLx,
0121 & 1-OLy:sNy+OLy,nSx,nSy)
0122 #endif /* ALLOW_DIAGNOSTICS */
0123
0124 _RL epsilon_H
0125 _RL stic_addMass(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
0126
0127 #ifdef ALLOW_GENTIM2D_CONTROL
0128 INTEGER iarr
0129 _RL xx_shifwflx_loc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
0130 #endif
0131
0132 #ifdef ALLOW_AUTODIFF_TAMC
0133
0134
0135
0136 INTEGER tkey, ikey, kkey
0137 #endif
0138
0139
0140
0141 #ifdef ALLOW_AUTODIFF_TAMC
0142
0143 #endif
0144
0145
0146
0147 epsilon_H = 1. _d -03
0148
0149
0150 thermalConductionDistance = 100.0 _d 0
0151 thermalConductionTemp = -20.0 _d 0
0152 icfgridareaFrac = 1.0 _d 0
0153 insituT = 0. _d 0
0154
0155
0156
0157
0158 DO bj = myByLo(myThid), myByHi(myThid)
0159 DO bi = myBxLo(myThid), myBxHi(myThid)
0160 DO j = 1-OLy,sNy+OLy
0161 DO i = 1-OLx,sNx+OLx
0162 shelfIceHeatFlux (i,j,bi,bj) = 0. _d 0
0163 shelfIceFreshWaterFlux (i,j,bi,bj) = 0. _d 0
0164 sticfHeatFlux (i,j,bi,bj) = 0. _d 0
0165 sticfFreshWaterFlux (i,j,bi,bj) = 0. _d 0
0166 #ifdef ALLOW_GENTIM2D_CONTROL
0167 xx_shifwflx_loc (i,j,bi,bj) = 0. _d 0
0168 #endif /* ALLOW_GENTIM2D_CONTROL */
0169 #ifndef ALLOW_SHITRANSCOEFF_3D
0170 shiTransCoeffS(i,j,bi,bj) = SHELFICEsaltToHeatRatio *
0171 & shiTransCoeffT(i,j,bi,bj)
0172 #endif
0173 ENDDO
0174 ENDDO
0175 DO k = 1, Nr
0176 DO j = 1-OLy,sNy+OLy
0177 DO i = 1-OLx,sNx+OLx
0178 #ifdef ALLOW_SHITRANSCOEFF_3D
0179 shiTransCoeffS3d(i,j,k,bi,bj) = SHELFICEsaltToHeatRatio *
0180 & shiTransCoeffT3d(i,j,k,bi,bj)
0181 #endif
0182 icfHeatFlux(i,j,k,bi,bj) = 0. _d 0
0183 icfFreshWaterFlux(i,j,k,bi,bj) = 0. _d 0
0184 stic_gT(i,j,k,bi,bj) = 0. _d 0
0185 stic_gS(i,j,k,bi,bj) = 0. _d 0
0186 stic_addMass(i,j,k,bi,bj) = 0. _d 0
0187 ENDDO
0188 ENDDO
0189 ENDDO
0190 #ifdef ALLOW_DIAGNOSTICS
0191 IF ( useDiagnostics ) THEN
0192 DO j = 1-OLy,sNy+OLy
0193 DO i = 1-OLx,sNx+OLx
0194 tmpDiagShelficeForcingT(i,j,bi,bj) = 0. _d 0
0195 tmpDiagShelficeForcingS(i,j,bi,bj) = 0. _d 0
0196 ENDDO
0197 ENDDO
0198 DO k = 1, Nr
0199 DO j = 1-OLy,sNy+OLy
0200 DO i = 1-OLx,sNx+OLx
0201 tmpDiagIcfForcingT(i,j,k,bi,bj) = 0. _d 0
0202 tmpDiagIcfForcingS(i,j,k,bi,bj) = 0. _d 0
0203 #ifdef NONLIN_FRSURF
0204 tmpDiag (i,j,k,bi,bj) = 0. _d 0
0205 tmpDiagT (i,j,k,bi,bj) = 0. _d 0
0206 #endif
0207 ENDDO
0208 ENDDO
0209 ENDDO
0210 ENDIF
0211 #endif /* ALLOW_DIAGNOSTICS */
0212 ENDDO
0213 ENDDO
0214
0215 #if (defined ALLOW_CTRL && defined ALLOW_GENTIM2D_CONTROL)
0216 IF ( useCTRL ) THEN
0217 DO iarr = 1, maxCtrlTim2D
0218 IF (xx_gentim2d_file(iarr)(1:11).EQ.'xx_shifwflx') THEN
0219 DO bj = myByLo(myThid),myByHi(myThid)
0220 DO bi = myBxLo(myThid),myBxHi(myThid)
0221 DO j = 1,sNy
0222 DO i = 1,sNx
0223 xx_shifwflx_loc(i,j,bi,bj)=xx_gentim2d(i,j,bi,bj,iarr)
0224 ENDDO
0225 ENDDO
0226 ENDDO
0227 ENDDO
0228 ENDIF
0229 ENDDO
0230 ENDIF
0231 #endif
0232
0233 DO bj = myByLo(myThid), myByHi(myThid)
0234 DO bi = myBxLo(myThid), myBxHi(myThid)
0235 #ifdef ALLOW_AUTODIFF_TAMC
0236 tkey = bi + (bj-1)*nSx + (ikey_dynamics-1)*nSx*nSy
0237 #endif
0238
0239
0240 DO j = 1-OLy+1,sNy+OLy-1
0241 DO i = 1-OLx+1,sNx+OLx-1
0242
0243 #ifdef ALLOW_AUTODIFF_TAMC
0244 ikey = i + OLx + (j+OLy-1)*(sNx+2*OLx)
0245 & + (tkey-1)*(sNx+2*OLx)*(sNy+2*OLy)
0246 #endif
0247
0248
0249 FRONT_K = kIcf(i,j,bi,bj)
0250
0251
0252 IF (FRONT_K .GT. 0) THEN
0253
0254
0255 DO k = 1, FRONT_K
0256
0257
0258
0259
0260
0261 DO SI = 1,4
0262 CURI=CURI_ARR(i,j,bi,bj,SI)
0263 CURJ=CURJ_ARR(i,j,bi,bj,SI)
0264 iceFrontWidth=sticfWidth_arr(i,j,bi,bj,SI)
0265
0266
0267
0268
0269 iceFrontCellThickness = 0. _d 0
0270 IF(iceFrontWidth.NE.0. _d 0)
0271 & iceFrontCellThickness = RA(CURI,CURJ,bi,bj)
0272 & /iceFrontWidth
0273 iceFrontFaceArea = drF(k)*iceFrontWidth
0274
0275
0276 IF (_hFacC(CURI,CURJ,k,bi,bj) .GT. zeroRL) THEN
0277
0278
0279
0280
0281
0282
0283
0284
0285
0286
0287
0288
0289
0290
0291
0292 ice_bottom_Z_C = max(rF(k+1),
0293 & min(Ro_surf(i,j, bi,bj), rF(k)))
0294
0295
0296
0297
0298
0299
0300
0301
0302
0303
0304 wet_top_Z_N = max(rF(k+1),
0305 & min(Ro_surf(CURI,CURJ, bi,bj), rF(k)))
0306
0307
0308
0309
0310
0311
0312
0313
0314
0315
0316
0317
0318
0319
0320
0321 wet_bottom_Z_N = min(rF(k),
0322 & max(R_low(CURI,CURJ, bi,bj), rF(k+1)))
0323
0324
0325
0326
0327
0328
0329
0330
0331
0332
0333
0334
0335 IF (ice_bottom_Z_C .GT. wet_bottom_Z_N) THEN
0336 iceFrontWetContact_Z_max = ice_bottom_Z_C
0337 ELSE
0338 iceFrontWetContact_Z_max = wet_bottom_Z_N
0339 ENDIF
0340
0341
0342
0343
0344
0345
0346
0347
0348
0349
0350 iceFrontVertContactFrac =
0351 & (wet_top_Z_N - iceFrontWetContact_Z_max)/ drF(k)
0352
0353
0354
0355
0356 IF (iceFrontVertContactFrac .GT. epsilon_H) THEN
0357 tLoc = theta(CURI,CURJ,k,bi,bj)
0358 sLoc = MAX(salt(CURI,CURJ,k,bi,bj), zeroRL)
0359
0360
0361
0362
0363 pLoc = 0.5 _d 0 * ABS(wet_top_Z_N +
0364 & iceFrontWetContact_Z_max)
0365 #ifdef ALLOW_SHITRANSCOEFF_3D
0366 gammaTLoc = shiTransCoeffT3d(CURI,CURJ,k,bi,bj)
0367 gammaSLoc = shiTransCoeffS3d(CURI,CURJ,k,bi,bj)
0368 #else
0369 gammaTLoc = shiTransCoeffT(CURI,CURJ,bi,bj)
0370 gammaSLoc = shiTransCoeffS(CURI,CURJ,bi,bj)
0371 #endif
0372
0373 CALL STIC_SOLVE4FLUXES(
0374 I tLoc, sLoc, pLoc,
0375 I gammaTLoc, gammaSLoc,
0376 I thermalConductionDistance,
0377 I thermalConductionTemp,
0378 O tmpHeatFlux, tmpFWFLX,
0379 O tmpForcingT, tmpForcingS,
0380 O insituT,
0381 I bi, bj, myTime, myIter, myThid )
0382 #ifdef ALLOW_AUTODIFF_TAMC
0383 kkey = k + (ikey-1)*Nr
0384
0385 #endif
0386
0387
0388
0389
0390
0391
0392
0393
0394
0395
0396
0397
0398
0399
0400 icfgridareaFrac =
0401 & iceFrontFaceArea/RA(CURI,CURJ,bi,bj)
0402
0403
0404 tmpHeatFluxscaled =
0405 & tmpHeatFlux*iceFrontVertContactFrac
0406 & *icfgridareaFrac
0407 icfHeatFlux(CURI,CURJ,k,bi,bj) =
0408 & icfHeatFlux(CURI,CURJ,k,bi,bj) +
0409 & tmpHeatFluxscaled
0410
0411
0412 tmpFWFLXscaled =
0413 & tmpFWFLX*iceFrontVertContactFrac
0414 & *icfgridareaFrac
0415 icfFreshWaterFlux(CURI,CURJ,k,bi,bj) =
0416 & icfFreshWaterFlux(CURI,CURJ,k,bi,bj) +
0417 & tmpFWFLXscaled
0418
0419
0420
0421
0422
0423
0424
0425 if(k.GE.kTopC(i,j,bi,bj))then
0426 if(RA(CURI,CURJ,bi,bj).NE.0. _d 0)then
0427 sticfHeatFlux(CURI,CURJ,bi,bj) =
0428 & sticfHeatFlux(CURI,CURJ,bi,bj) +
0429 & tmpHeatFluxscaled
0430 sticfFreshWaterFlux(CURI,CURJ,bi,bj) =
0431 & sticfFreshWaterFlux(CURI,CURJ,bi,bj) +
0432 & tmpFWFLXscaled
0433 endif
0434 endif
0435
0436
0437
0438
0439
0440
0441
0442 IF (iceFrontCellThickness .NE. 0. _d 0) THEN
0443
0444 stic_gT(CURI,CURJ,k,bi,bj) =
0445 & stic_gT(CURI,CURJ,k,bi,bj) +
0446 & tmpForcingT/iceFrontCellThickness*
0447 & iceFrontVertContactFrac*
0448 & _recip_hFacC(CURI,CURJ,k,bi,bj)
0449
0450 stic_gS(CURI,CURJ,k,bi,bj) =
0451 & stic_gS(CURI,CURJ,k,bi,bj) +
0452 & tmpForcingS/iceFrontCellThickness*
0453 & iceFrontVertContactFrac*
0454 & _recip_hFacC(CURI,CURJ,k,bi,bj)
54f5c8e380 Jean*0455 #ifdef ALLOW_DIAGNOSTICS
0456 IF ( useDiagnostics ) THEN
0457 tmpDiagIcfForcingT(CURI,CURJ,k,bi,bj) =
0458 & tmpDiagIcfForcingT(CURI,CURJ,k,bi,bj) +
0459 & tmpForcingT/iceFrontCellThickness*
0460 & iceFrontVertContactFrac*
0461 & drF(k)
0462 tmpDiagIcfForcingS(CURI,CURJ,k,bi,bj) =
0463 & tmpDiagIcfForcingS(CURI,CURJ,k,bi,bj) +
0464 & tmpForcingS/iceFrontCellThickness*
0465 & iceFrontVertContactFrac*
0466 & drF(k)
0467 ENDIF
0468 #endif /* ALLOW_DIAGNOSTICS */
00f81e6785 Ou W*0469 ENDIF /* iceFrontCellThickness */
0470
0471 stic_addMass(CURI,CURJ,k,bi,bj) =
0472 & stic_addMass(CURI,CURJ,k,bi,bj) -
0473 & tmpFWFLX*iceFrontFaceArea*
0474 & iceFrontVertContactFrac
0475 #ifdef ALLOW_DIAGNOSTICS
0476 IF ( useDiagnostics ) THEN
0477 #ifdef NONLIN_FRSURF
0478 IF ( (nonlinFreeSurf.GT.0 .OR. usingPCoords)
0479 & .AND. useRealFreshWaterFlux ) THEN
0480 tmpDiagT(CURI,CURJ,k,bi,bj) =
0481 & stic_addMass(CURI,CURJ,k, bi,bj)/
0482 & RA(CURI,CURJ,bi,bj)*
0483 & insituT*HeatCapacity_Cp
0484 ENDIF
0485 #endif /* NONLIN_FRSURF */
0486 ENDIF
0487 #endif /* ALLOW_DIAGNOSTICS */
0488 ENDIF /* iceFrontVertContactFrac */
0489 ENDIF /* hFacC(CURI,CURJ,k,bi,bj) */
0490 ENDDO /* SI loop for adjacent cells */
0491 ENDDO /* k Loop */
0492 ENDIF /* FRONT_K */
0493
0494
0495 k = kTopC(i,j,bi,bj)
0496
0497
0498
0499
0500 IF (k .GT. 0) THEN
0501
0502 pLoc = 0. _d 0
0503 tLoc = 0. _d 0
0504 sLoc = 0. _d 0
0505 gammaTLoc = 0. _d 0
0506 gammaSLoc = 0. _d 0
0507
0508
0509
0510
0511
0512
0513 pLoc = ABS(R_shelfIce(i,j,bi,bj))
0514 tLoc = theta(i,j,k,bi,bj)
0515 sLoc = MAX(salt(i,j,k,bi,bj), zeroRL)
0516 #ifdef ALLOW_SHITRANSCOEFF_3D
0517 gammaTLoc = shiTransCoeffT3d(i,j,k,bi,bj)
0518 gammaSLoc = shiTransCoeffS3d(i,j,k,bi,bj)
0519 #else
0520 gammaTLoc = shiTransCoeffT(i,j,bi,bj)
0521 gammaSLoc = shiTransCoeffS(i,j,bi,bj)
0522 #endif
0523 CALL STIC_SOLVE4FLUXES(
0524 I tLoc, sLoc, pLoc,
0525 I gammaTLoc, gammaSLoc,
0526 I pLoc, thermalConductionTemp,
0527 O tmpHeatFlux, tmpFWFLX,
0528 O tmpForcingT, tmpForcingS,
0529 O insituT,
0530 I bi, bj, myTime, myIter, myThid )
0531
0532 shelficeHeatFlux(i,j,bi,bj) = tmpHeatFlux
0533 #ifdef ALLOW_GENTIM2D_CONTROL
0534 & - xx_shifwflx_loc(i,j,bi,bj)*SHELFICElatentHeat
0535 #endif /* ALLOW_GENTIM2D_CONTROL */
0536
0537 shelfIceFreshWaterFlux(i,j,bi,bj) = tmpFWFLX
0538 #ifdef ALLOW_GENTIM2D_CONTROL
0539 & + xx_shifwflx_loc(i,j,bi,bj)
0540 #endif /* ALLOW_GENTIM2D_CONTROL */
0541
0542
0543 sticfHeatFlux(i,j,bi,bj) =
0544 & sticfHeatFlux(i,j,bi,bj) +
0545 & shelficeHeatFlux(i,j,bi,bj)
0546 sticfFreshWaterFlux(i,j,bi,bj) =
0547 & sticfFreshWaterFlux(i,j,bi,bj) +
0548 & shelfIceFreshWaterFlux(i,j,bi,bj)
0549
0550 stic_gT(i,j,k,bi,bj) = stic_gT(i,j,k,bi,bj) +
0551 & tmpForcingT*recip_drF(k)* _recip_hFacC(i,j,k,bi,bj)
0552
0553 stic_gS(i,j,k,bi,bj) = stic_gS(i,j,k,bi,bj) +
0554 & tmpForcingS*recip_drF(k)* _recip_hFacC(i,j,k,bi,bj)
0555
0556 stic_addMass(i,j,k,bi,bj) = stic_addMass(i,j,k,bi,bj) -
0557 & tmpFWFLX*RA(i,j,bi,bj)
0558 #ifdef ALLOW_DIAGNOSTICS
0559 IF ( useDiagnostics ) THEN
54f5c8e380 Jean*0560 tmpDiagShelficeForcingT(i,j,bi,bj) = tmpForcingT
0561 tmpDiagShelficeForcingS(i,j,bi,bj) = tmpForcingS
00f81e6785 Ou W*0562 #ifdef NONLIN_FRSURF
0563 IF ( (nonlinFreeSurf.GT.0 .OR. usingPCoords)
0564 & .AND. useRealFreshWaterFlux ) THEN
0565 tmpDiagT(i,j,k,bi,bj) =
0566 & stic_addMass(i,j,k,bi,bj) / rA(i,j,bi,bj)*
0567 & insituT*HeatCapacity_Cp
0568 ENDIF
0569 #endif /* NONLIN_FRSURF */
0570 ENDIF
0571 #endif /* ALLOW_DIAGNOSTICS */
0572 ENDIF /* Shelf k > 0 */
0573 ENDDO /* i */
0574 ENDDO /* j */
0575 ENDDO /* bi */
0576 ENDDO /* bj */
0577
0578 DO bj = myByLo(myThid), myByHi(myThid)
0579 DO bi = myBxLo(myThid), myBxHi(myThid)
0580 DO j = 1-OLy+1,sNy+OLy-1
0581 DO i = 1-OLx+1,sNx+OLx-1
0582 DO k = 1,Nr
0583 addMass(i,j,k,bi,bj) = addMass(i,j,k,bi,bj) +
0584 & stic_addMass(i,j,k,bi,bj)
0585 ENDDO /* k */
0586 ENDDO /* i */
0587 ENDDO /* j */
0588 ENDDO /* bi */
0589 ENDDO /* bj */
0590
0591 #ifdef ALLOW_ADDFLUID
0592 IF ( selectAddFluid.NE.0 .AND. salt_addMass.NE.UNSET_RL ) THEN
0593 IF ( ( selectAddFluid.GE.1 .AND. nonlinFreeSurf.GT.0 )
0594 & .OR. convertFW2Salt.EQ.-1. _d 0 ) THEN
0595 DO bj = myByLo(myThid), myByHi(myThid)
0596 DO bi = myBxLo(myThid), myBxHi(myThid)
0597 DO j=0,sNy+1
0598 DO i=0,sNx+1
0599 DO k=1,Nr
0600 stic_gT(i,j,k,bi,bj) = stic_gT(i,j,k,bi,bj)
0601 & - stic_addMass(i,j,k,bi,bj)*mass2rUnit
0602 & *( temp_addMass - theta(i,j,k,bi,bj) )
0603 & *recip_rA(i,j,bi,bj)
0604 & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
0605
0606 stic_gS(i,j,k,bi,bj) = stic_gS(i,j,k,bi,bj)
0607 & - stic_addMass(i,j,k,bi,bj)*mass2rUnit
0608 & *( salt_addMass - salt(i,j,k,bi,bj) )
0609 & *recip_rA(i,j,bi,bj)
0610 & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
0611
0612 ENDDO /* k */
0613 ENDDO /* i */
0614 ENDDO /* j */
0615 ENDDO /* bi */
0616 ENDDO /* bj */
0617 ELSE
0618 DO bj = myByLo(myThid), myByHi(myThid)
0619 DO bi = myBxLo(myThid), myBxHi(myThid)
0620 DO j=0,sNy+1
0621 DO i=0,sNx+1
0622 DO k=1,Nr
0623 stic_gT(i,j,k,bi,bj) = stic_gT(i,j,k,bi,bj)
0624 & - stic_addMass(i,j,k,bi,bj)*mass2rUnit
0625 & *( temp_addMass - tRef(k) )
0626 & *recip_rA(i,j,bi,bj)
0627 & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
0628
0629 stic_gS(i,j,k,bi,bj) = stic_gS(i,j,k,bi,bj)
0630 & - stic_addMass(i,j,k,bi,bj)*mass2rUnit
0631 & *( salt_addMass - sRef(k) )
0632 & *recip_rA(i,j,bi,bj)
0633 & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
0634
0635 ENDDO /* k */
0636 ENDDO /* i */
0637 ENDDO /* j */
0638 ENDDO /* bi */
0639 ENDDO /* bj */
0640 ENDIF
0641 ENDIF
0642 #endif /* ALLOW_ADDFLUID */
0643
0644
0645 #ifndef ALLOW_AUTODIFF
0646
0647 DO bj = myByLo(myThid), myByHi(myThid)
0648 DO bi = myBxLo(myThid), myBxHi(myThid)
0649 DO j = 1-OLy, sNy+OLy
0650 DO i = 1-OLx, sNx+OLx
0651 shelficeLoadAnomaly(i,j,bi,bj) = gravity
0652 & *( shelficeMass(i,j,bi,bj) + rhoConst*Ro_surf(i,j,bi,bj) )
0653 ENDDO
0654 ENDDO
0655 ENDDO
0656 ENDDO
0657
0658 #endif /* ndef ALLOW_AUTODIFF */
0659
0660 #ifdef ALLOW_DIAGNOSTICS
0661 IF ( useDiagnostics ) THEN
54f5c8e380 Jean*0662
00f81e6785 Ou W*0663 #ifdef NONLIN_FRSURF
0664 IF ( (nonlinFreeSurf.GT.0 .OR. usingPCoords)
0665 & .AND. useRealFreshWaterFlux ) THEN
0666 DO bj = myByLo(myThid), myByHi(myThid)
0667 DO bi = myBxLo(myThid), myBxHi(myThid)
0668 DO k = 1, Nr
0669 DO j=1,sNy
0670 DO i=1,sNx
0671 tmpDiag(i,j,k,bi,bj) =
0672 & stic_addMass(i,j,k,bi,bj)/RA(i,j,bi,bj)*
0673 & theta(i,j,k,bi,bj)*HeatCapacity_Cp
0674 ENDDO
0675 ENDDO
0676 ENDDO
0677 ENDDO /* bi */
0678 ENDDO /* bj */
0679 CALL DIAGNOSTICS_FILL( tmpDiag,'TFLUXMLT',0,Nr,0,1,1,myThid )
0680
0681 DO bj = myByLo(myThid), myByHi(myThid)
0682 DO bi = myBxLo(myThid), myBxHi(myThid)
0683 DO k = 1, Nr
0684 DO j=1,sNy
0685 DO i=1,sNx
0686 tmpDiag(i,j,k,bi,bj) =
0687 & stic_addMass(i,j,k,bi,bj)/RA(i,j,bi,bj)*salt(i,j,k,bi,bj)
0688 ENDDO
0689 ENDDO
0690 ENDDO
0691 ENDDO /* bi */
0692 ENDDO /* bj */
0693 CALL DIAGNOSTICS_FILL( tmpDiag,'SFLUXMLT',0,Nr,0,1,1,myThid )
0694
0695 CALL DIAGNOSTICS_FILL( tmpDiagT,'TFLUXMTI',0,Nr,0,1,1,myThid )
0696 ENDIF
0697 #endif /* NONLIN_FRSURF */
0698
0699 CALL DIAGNOSTICS_FILL_RS(shelfIceFreshWaterFlux,'SHIfwFlx',
0700 & 0,1,0,1,1,myThid)
0701 CALL DIAGNOSTICS_FILL_RS(shelfIceHeatFlux, 'SHIhtFlx',
0702 & 0,1,0,1,1,myThid)
0703
0704 #ifdef ALLOW_STEEP_ICECAVITY
54f5c8e380 Jean*0705 CALL DIAGNOSTICS_FILL( sticfFreshWaterFlux, 'SHIICFfw',
00f81e6785 Ou W*0706 & 0,1,0,1,1,myThid)
54f5c8e380 Jean*0707 CALL DIAGNOSTICS_FILL( sticfHeatFlux, 'SHIICFht',
00f81e6785 Ou W*0708 & 0,1,0,1,1,myThid)
0709 #endif
0710
0711 CALL DIAGNOSTICS_FILL(icfFreshWaterFlux, 'STCfwFlx',
0712 & 0,Nr,0,1,1,myThid)
0713 CALL DIAGNOSTICS_FILL(icfHeatFlux, 'STChtFlx',
0714 & 0,Nr,0,1,1,myThid)
0715
0716
0717 tmpFac = HeatCapacity_Cp*rUnit2mass
0718 CALL DIAGNOSTICS_SCALE_FILL(tmpDiagShelficeForcingT,tmpFac,1,
0719 & 'SHIForcT',0,1,0,1,1,myThid)
0720
0721 tmpFac = rUnit2mass
0722 CALL DIAGNOSTICS_SCALE_FILL(tmpDiagShelficeForcingS,tmpFac,1,
0723 & 'SHIForcS',0,1,0,1,1,myThid)
0724
0725
0726 tmpFac = HeatCapacity_Cp*rUnit2mass
0727 CALL DIAGNOSTICS_SCALE_FILL(tmpDiagIcfForcingT,tmpFac,1,
0728 & 'STCForcT',0,Nr,0,1,1,myThid)
0729
0730 tmpFac = rUnit2mass
0731 CALL DIAGNOSTICS_SCALE_FILL(tmpDiagIcfForcingS,tmpFac,1,
0732 & 'STCForcS',0,Nr,0,1,1,myThid)
0733
0734
0735 #ifdef ALLOW_SHITRANSCOEFF_3D
0736 CALL DIAGNOSTICS_FILL(shiTransCoeffT3d,'SHIgam3T',
0737 & 0,Nr,0,1,1,myThid)
0738 CALL DIAGNOSTICS_FILL(shiTransCoeffS3d,'SHIgam3S',
0739 & 0,Nr,0,1,1,myThid)
0740 #else
0741 CALL DIAGNOSTICS_FILL(shiTransCoeffT,'SHIgammT',
0742 & 0,1,0,1,1,myThid)
0743 CALL DIAGNOSTICS_FILL(shiTransCoeffS,'SHIgammS',
0744 & 0,1,0,1,1,myThid)
0745 #endif
0746
0747 #ifdef SHI_ALLOW_GAMMAFRICT
0748 IF ( SHELFICEuseGammaFrict )
0749 & CALL DIAGNOSTICS_FILL(uStarDiag,'SHIuStar',0,1,0,1,1,myThid)
0750 #endif /* SHI_ALLOW_GAMMAFRICT */
0751
0752 ENDIF
54f5c8e380 Jean*0753 #endif /* ALLOW_DIAGNOSTICS */
00f81e6785 Ou W*0754
0755 RETURN
0756 END