/[MITgcm]/MITgcm_contrib/shelfice_remeshing/MANUAL/code/shelfice_thermodynamics.F
ViewVC logotype

Diff of /MITgcm_contrib/shelfice_remeshing/MANUAL/code/shelfice_thermodynamics.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.2 by dgoldberg, Mon Sep 28 12:59:59 2015 UTC revision 1.3 by dgoldberg, Tue Oct 13 15:54:41 2015 UTC
# Line 82  C     a0, a1, a2, b, c0 Line 82  C     a0, a1, a2, b, c0
82  C     eps1, eps2, eps3, eps3a, eps4, eps5, eps6, eps7, eps8  C     eps1, eps2, eps3, eps3a, eps4, eps5, eps6, eps7, eps8
83  C     aqe, bqe, cqe, discrim, recip_aqe  C     aqe, bqe, cqe, discrim, recip_aqe
84  C     drKp1, recip_drLoc  C     drKp1, recip_drLoc
85        INTEGER I,J,K,Kp1        INTEGER I,J,K,Kp1,kp2
86        INTEGER bi,bj        INTEGER bi,bj
87        _RL tLoc(1:sNx,1:sNy)        _RL tLoc(1:sNx,1:sNy)
88        _RL sLoc(1:sNx,1:sNy)        _RL sLoc(1:sNx,1:sNy)
# Line 95  C     drKp1, recip_drLoc Line 95  C     drKp1, recip_drLoc
95        _RL freshWaterFlux, convertFW2SaltLoc        _RL freshWaterFlux, convertFW2SaltLoc
96        _RL a0, a1, a2, b, c0        _RL a0, a1, a2, b, c0
97        _RL eps1, eps2, eps3, eps3a, eps4, eps5, eps6, eps7, eps8        _RL eps1, eps2, eps3, eps3a, eps4, eps5, eps6, eps7, eps8
98        _RL cFac, rFac, dFac, fwflxFac        _RL cFac, rFac, dFac, fwflxFac, realfwFac
99        _RL aqe, bqe, cqe, discrim, recip_aqe        _RL aqe, bqe, cqe, discrim, recip_aqe
100        _RL drKp1, recip_drLoc        _RL drKp1, drKp2, recip_drLoc
101        _RL recip_latentHeat        _RL recip_latentHeat
102        _RL tmpFac        _RL tmpFac
103    
# Line 173  C     are we doing the conservative form Line 173  C     are we doing the conservative form
173        recip_Cp = 1. _d 0 / HeatCapacity_Cp        recip_Cp = 1. _d 0 / HeatCapacity_Cp
174        cFac = 0. _d 0        cFac = 0. _d 0
175        IF ( SHELFICEconserve ) cFac = 1. _d 0        IF ( SHELFICEconserve ) cFac = 1. _d 0
176    
177          realFWfac = 0. _d 0
178          IF ( SHELFICErealFWflux ) realFWfac = 1. _d 0
179  C     with "real fresh water flux" (affecting ETAN),  C     with "real fresh water flux" (affecting ETAN),
180  C     there is more to modify  C     there is more to modify
181        rFac = 1. _d 0        rFac = 1. _d 0
# Line 274  C--   underneath the ice Line 277  C--   underneath the ice
277  c         pLoc(I,J) = shelficeMass(I,J,bi,bj)*gravity*1. _d -4  c         pLoc(I,J) = shelficeMass(I,J,bi,bj)*gravity*1. _d -4
278            tLoc(I,J) = theta(I,J,K,bi,bj)            tLoc(I,J) = theta(I,J,K,bi,bj)
279            sLoc(I,J) = MAX(salt(I,J,K,bi,bj), zeroRL)            sLoc(I,J) = MAX(salt(I,J,K,bi,bj), zeroRL)
 c              print *, 'JJHA' , hFacC(I,J,K,bi,bj)            
280            IF ( .not.SHELFICEBoundaryLayer ) THEN            IF ( .not.SHELFICEBoundaryLayer ) THEN
281             uLoc(I,J) = recip_hFacC(I,J,K,bi,bj) *             uLoc(I,J) = recip_hFacC(I,J,K,bi,bj) *
282       &         ( uVel(I,  J,K,bi,bj) * _hFacW(I,  J,K,bi,bj)       &         ( uVel(I,  J,K,bi,bj) * _hFacW(I,  J,K,bi,bj)
# Line 283  c              print *, 'JJHA' , hFacC(I Line 285  c              print *, 'JJHA' , hFacC(I
285       &         ( vVel(I,  J,K,bi,bj) * _hFacS(I,  J,K,bi,bj)       &         ( vVel(I,  J,K,bi,bj) * _hFacS(I,  J,K,bi,bj)
286       &         + vVel(I,J+1,K,bi,bj) * _hFacS(I,J+1,K,bi,bj) )       &         + vVel(I,J+1,K,bi,bj) * _hFacS(I,J+1,K,bi,bj) )
287            ENDIF            ENDIF
 !           tmpfac =  
 !      &         ( _hFacS(I,J,  K,bi,bj) +  
 !      &           _hFacS(I,J+1,K,bi,bj) ) / 2.0  
 !           vLoc(I,J) =  
 !      &         ( vVel(I,J,  K,bi,bj) * _hFacS(I,J,  K,bi,bj)  
 !      &         + vVel(I,J+1,K,bi,bj) * _hFacS(I,J+1,K,bi,bj) )  
 !           if (tmpfac .gt. 0.0) then  
 !                if (j.lt.100) then  
 !                print *, 1./tmpfac, recip_hFacC(I,J,K,bi,bj)  
 !                endif  
 !                vLoc(I,J) = vLoc(I,J) / tmpfac  
 !           else  
 !                vLoc(I,J) = 0.0  
 !           endif  
288           ENDDO           ENDDO
289          ENDDO          ENDDO
290    
291    !         IF ( SHELFICEBoundaryLayer ) THEN
292    !          DO J = 1, sNy+1
293    !           DO I = 1, sNx+1
294    !
295    !            K = ksurfW(I,J,bi,bj)
296    !            Kp1 = K+1
297    !            Kp2 = K+2
298    !
299    !            IF (ShelficeThickBoundaryLayer .and.
300    !      &      (K.ne.0.and.K.LT.Nr-1)) THEN
301    !
302    !             drKp1 = drF(K)*( 1.5 - _hFacW(I,J,K,bi,bj) )
303    !             drKp2 = drKp1 - drF(kp1)*_hFacW(I,J,kp1,bi,bj)
304    !             drKp2 = MAX( drKp2, 0. _d 0)
305    !             drKp2 = MIN( drKp2,
306    !      &        drF(kp2)*_hFacW(I,J,kp2,bi,bj))
307    !             drKp1 = drKp1 - drKp2
308    !             drKp1 = MAX( drKp1, 0. _d 0)
309    !             recip_drLoc = 1. _d 0 /
310    !      &           (drF(K)*_hFacW(I,J,K,bi,bj)+drKp1+drKp2)            
311    !             u_topdr(I,J,bi,bj) =
312    !      &       (drF(K)*_hFacW(I,J,K,bi,bj)*uVel(I,J,K,bi,bj) +
313    !      &        drKp1*uVel(I,J,Kp1,bi,bj)) * recip_drLoc
314    !             u_topdr(I,J,bi,bj) = u_topdr(I,J,bi,bj) +
315    !      &        drKp2 * uVel(I,J,Kp2,bi,bj) * recip_drLoc
316    !
317    !            ELSEIF ( (K .NE. 0 .AND. K.EQ.Nr-1) .OR.
318    !      &      (.not.SHELFICEthickboundarylayer.AND.
319    !      &       (K .NE. 0 .AND. K .LT. Nr) ) ) THEN
320    !
321    !             drKp1 = drF(K)*(1. _d 0-_hFacW(I,J,K,bi,bj))
322    !             drKp1 = max (drKp1, 0. _d 0)
323    !             recip_drLoc = 1.0 /
324    !      &       (drF(K)*_hFacW(I,J,K,bi,bj)+drKp1)
325    !             u_topdr(I,J,bi,bj) =
326    !      &       (drF(K)*_hFacW(I,J,K,bi,bj)*uVel(I,J,K,bi,bj) +
327    !      &        drKp1*uVel(I,J,Kp1,bi,bj))
328    !      &      * recip_drLoc
329    !
330    !            ELSE
331    !
332    !             u_topdr(I,J,bi,bj) = 0. _d 0
333    !
334    !            ENDIF
335    !            
336    !            K = ksurfS(I,J,bi,bj)
337    !            Kp1 = K+1
338    !            Kp2 = K+2
339    !
340    !            IF (ShelficeThickBoundaryLayer .and.
341    !      &      (K.ne.0.and.K.LT.Nr-1)) THEN
342    !
343    !             drKp1 = drF(K)*( 1.5 - _hFacS(I,J,K,bi,bj) )
344    !             drKp2 = drKp1 - drF(kp1)*_hFacS(I,J,kp1,bi,bj)
345    !             drKp2 = MAX( drKp2, 0. _d 0)
346    !             drKp2 = MIN( drKp2,
347    !      &        drF(kp2)*_hFacS(I,J,kp2,bi,bj))
348    !             drKp1 = drKp1 - drKp2
349    !             drKp1 = MAX( drKp1, 0. _d 0)
350    !             recip_drLoc = 1. _d 0 /
351    !      &           (drF(K)*_hFacS(I,J,K,bi,bj)+drKp1+drKp2)            
352    !             v_topdr(I,J,bi,bj) =
353    !      &       (drF(K)*_hFacS(I,J,K,bi,bj)*vVel(I,J,K,bi,bj) +
354    !      &        drKp1*vVel(I,J,Kp1,bi,bj)) * recip_drLoc
355    !             v_topdr(I,J,bi,bj) = v_topdr(I,J,bi,bj) +
356    !      &        drKp2 * vVel(I,J,Kp2,bi,bj) * recip_drLoc
357    !
358    !            ELSEIF ( (K .NE. 0 .AND. K.EQ.Nr-1) .OR.
359    !      &      ((.NOT.SHELFICEthickboundarylayer).AND.
360    !      &       (K .NE. 0 .AND. K .LT. Nr) ) ) THEN
361    !
362    !             drKp1 = drF(K)*(1. _d 0-_hFacS(I,J,K,bi,bj))
363    !             drKp1 = max (drKp1, 0. _d 0)
364    !             recip_drLoc = 1.0 /
365    !      &       (drF(K)*_hFacS(I,J,K,bi,bj)+drKp1)
366    !             v_topdr(I,J,bi,bj) =
367    !      &       (drF(K)*_hFacS(I,J,K,bi,bj)*vVel(I,J,K,bi,bj) +
368    !      &        drKp1*vVel(I,J,Kp1,bi,bj))
369    !      &      * recip_drLoc
370    !
371    !            ELSE
372    !
373    !             v_topdr(I,J,bi,bj) = 0. _d 0
374    !
375    !            ENDIF
376    !
377    !           ENDDO
378    !          ENDDO
379    !         ENDIF
380    
381          IF ( SHELFICEBoundaryLayer ) THEN          IF ( SHELFICEBoundaryLayer ) THEN
382           DO J = 1, sNy+1           DO J = 1, sNy+1
383            DO I = 1, sNx+1            DO I = 1, sNx+1
# Line 374  C--   lower cell may not be as thick as Line 452  C--   lower cell may not be as thick as
452           ENDDO           ENDDO
453          ENDIF          ENDIF
454    
455    
456          IF ( SHELFICEBoundaryLayer ) THEN          IF ( SHELFICEBoundaryLayer ) THEN
457           DO J = 1, sNy           DO J = 1, sNy
458            DO I = 1, sNx            DO I = 1, sNx
# Line 600  C--   compute surface tendencies Line 679  C--   compute surface tendencies
679       &           ( shiTransCoeffT(i,j,bi,bj)       &           ( shiTransCoeffT(i,j,bi,bj)
680       &           - cFac*shelfIceFreshWaterFlux(I,J,bi,bj)*mass2rUnit )       &           - cFac*shelfIceFreshWaterFlux(I,J,bi,bj)*mass2rUnit )
681       &           * ( thetaFreeze - tLoc(I,J) )       &           * ( thetaFreeze - tLoc(I,J) )
682         &           - realFWfac*shelfIceFreshWaterFlux(I,J,bi,bj)*
683         &             mass2rUnit*
684         &             ( tLoc(I,J) - theta(I,J,K,bi,bj) )
685              shelficeForcingS(i,j,bi,bj) =              shelficeForcingS(i,j,bi,bj) =
686       &           ( shiTransCoeffS(i,j,bi,bj)       &           ( shiTransCoeffS(i,j,bi,bj)
687       &           - cFac*shelfIceFreshWaterFlux(I,J,bi,bj)*mass2rUnit )       &           - cFac*shelfIceFreshWaterFlux(I,J,bi,bj)*mass2rUnit )
688       &           * ( saltFreeze - sLoc(I,J) )       &           * ( saltFreeze - sLoc(I,J) )
689         &           - realFWfac*shelfIceFreshWaterFlux(I,J,bi,bj)*
690         &             mass2rUnit*
691         &             ( sLoc(I,J) - salt(I,J,K,bi,bj) )
692             ELSE             ELSE
693              shelfIceHeatFlux      (I,J,bi,bj) = 0. _d 0              shelfIceHeatFlux      (I,J,bi,bj) = 0. _d 0
694              shelfIceFreshWaterFlux(I,J,bi,bj) = 0. _d 0              shelfIceFreshWaterFlux(I,J,bi,bj) = 0. _d 0
# Line 618  C     endif (not) useISOMIPTD Line 703  C     endif (not) useISOMIPTD
703        ENDDO        ENDDO
704    
705        IF (SHELFICEMassStepping) THEN        IF (SHELFICEMassStepping) THEN
706         CALL SHELFICE_STEP_ICEMASS( myTime, myIter, myThid )  !       CALL SHELFICE_STEP_ICEMASS( myTime, myIter, myThid )
707        ENDIF        ENDIF
708    
709  C--  Calculate new loading anomaly (in case the ice-shelf mass was updated)  C--  Calculate new loading anomaly (in case the ice-shelf mass was updated)
# Line 665  C     Friction velocity Line 750  C     Friction velocity
750  #endif /* ALLOW_DIAGNOSTICS */  #endif /* ALLOW_DIAGNOSTICS */
751    
752  #endif /* ALLOW_SHELFICE */  #endif /* ALLOW_SHELFICE */
         
753        RETURN        RETURN
   
754        END        END
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.3
  ViewVC Help
Powered by ViewVC 1.1.22