/[MITgcm]/MITgcm/pkg/seaice/seaice_growth.F
ViewVC logotype

Diff of /MITgcm/pkg/seaice/seaice_growth.F

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

revision 1.86 by jmc, Thu Oct 14 16:05:59 2010 UTC revision 1.87 by gforget, Thu Oct 14 16:29:44 2010 UTC
# Line 1  Line 1 
1    #define SEAICE_GROWTH_LEGACY
2  C $Header$  C $Header$
3  C $Name$  C $Name$
4    
# Line 88  c     prone to generate SNOW (FRWfromSNW Line 89  c     prone to generate SNOW (FRWfromSNW
89  c ICE/SNOW stocks tendencies associated with the various melt/freeze processes  c ICE/SNOW stocks tendencies associated with the various melt/freeze processes
90        _RL d_AREAbyATM         (1:sNx,1:sNy)        _RL d_AREAbyATM         (1:sNx,1:sNy)
91  c  c
92          _RL d_HEFFbyNEG         (1:sNx,1:sNy)
93        _RL d_HEFFbyOCNonICE    (1:sNx,1:sNy)        _RL d_HEFFbyOCNonICE    (1:sNx,1:sNy)
94        _RL d_HEFFbyATMonOCN    (1:sNx,1:sNy)        _RL d_HEFFbyATMonOCN    (1:sNx,1:sNy)
95        _RL d_HEFFbyFLOODING    (1:sNx,1:sNy)        _RL d_HEFFbyFLOODING    (1:sNx,1:sNy)
96  c  c
97          _RL d_HSNWbyNEG         (1:sNx,1:sNy)
98        _RL d_HSNWbyATMonSNW    (1:sNx,1:sNy)        _RL d_HSNWbyATMonSNW    (1:sNx,1:sNy)
99        _RL d_HSNWbyOCNonSNW    (1:sNx,1:sNy)        _RL d_HSNWbyOCNonSNW    (1:sNx,1:sNy)
100        _RL d_HSNWbyRAIN        (1:sNx,1:sNy)        _RL d_HSNWbyRAIN        (1:sNx,1:sNy)
# Line 102  C     actual ice thickness with upper an Line 105  C     actual ice thickness with upper an
105        _RL heffActual          (1:sNx,1:sNy)        _RL heffActual          (1:sNx,1:sNy)
106  C     actual snow thickness  C     actual snow thickness
107        _RL hsnowActual         (1:sNx,1:sNy)        _RL hsnowActual         (1:sNx,1:sNy)
108    
109    C     AREA_PRE :: hold sea-ice fraction field before any seaice-thermo update
110          _RL AREApreTH           (1:sNx,1:sNy)
111          _RL HEFFpreTH           (1:sNx,1:sNy)
112          _RL HSNWpreTH           (1:sNx,1:sNy)
113    
114  C     wind speed  C     wind speed
115        _RL UG                  (1:sNx,1:sNy)        _RL UG                  (1:sNx,1:sNy)
116        _RL SPEED_SQ        _RL SPEED_SQ
117  C temporary variables used to compute/regularize "actual" thicknesses  C temporary variables used to compute/regularize "actual" thicknesses
118        _RL areaSup,areaMin,hiceMin        _RL areaMin,hiceMin
119    
120  c temporary variables available for the various computations  c temporary variables available for the various computations
121        _RL tmpscal1, tmpscal2, tmpscal3, tmpscal4        _RL tmpscal1, tmpscal2, tmpscal3, tmpscal4
122        _RL tmparr1             (1:sNx,1:sNy)        _RL tmparr1             (1:sNx,1:sNy)
123    
 C auxillary variables used for specific processes  
       _RL snowEnergy  
124    
125  #ifdef ALLOW_SEAICE_FLOODING  #ifdef ALLOW_SEAICE_FLOODING
126        _RL hDraft        _RL hDraft
# Line 132  C auxillary variables used for specific Line 139  C auxillary variables used for specific
139        _RL a_QSWbyATMmult_cover(1:sNx,1:sNy)        _RL a_QSWbyATMmult_cover(1:sNx,1:sNy)
140  #endif  #endif
141    
 #ifdef SEAICE_AGE  
 C     old_AREA :: hold sea-ice fraction field before any seaice-thermo update  
       _RL old_AREA            (1:sNx,1:sNy)  
 # ifdef SEAICE_AGE_VOL  
 C     old_HEFF :: hold sea-ice effective thickness field before any seaice-thermo update  
       _RL old_HEFF            (1:sNx,1:sNy)  
       _RL age_actual  
 # endif /* SEAICE_AGE_VOL */  
 #endif /* SEAICE_AGE */  
   
142  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
143        _RL DIAGarray     (1:sNx,1:sNy)        _RL DIAGarray     (1:sNx,1:sNy)
144        LOGICAL  DIAGNOSTICS_IS_ON        LOGICAL  DIAGNOSTICS_IS_ON
# Line 260  c Line 257  c
257          ENDDO          ENDDO
258  #endif  #endif
259    
260  #ifdef SEAICE_AGE  
261  C     store the initial ice fraction over the domain  c ===================================================================
262    c =============PART 1: regularization, after advdiff, etc.===========
263    c ===================================================================
264    
265    #ifndef SEAICE_GROWTH_LEGACY
266    
267    c 1) treat the case of negative values:
268    c
269    #ifdef ALLOW_AUTODIFF_TAMC
270    CADJ STORE heff(:,:,bi,bj)  = comlev1_bibj,
271    CADJ &                        key = iicekey, byte = isbyte
272    CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,
273    CADJ &                        key = iicekey, byte = isbyte
274    CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
275    CADJ &                        key = iicekey, byte = isbyte
276    #endif /* ALLOW_AUTODIFF_TAMC */
277          DO J=1,sNy          DO J=1,sNy
278           DO I=1,sNx           DO I=1,sNx
279             old_AREA(i,j) = AREA(I,J,bi,bj)            d_HEFFbyNEG(I,J)=MAX(-HEFF(I,J,bi,bj),0._d 0)
280  # ifdef SEAICE_AGE_VOL            HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj)+d_HEFFbyNEG(I,J)
281             old_HEFF(i,j) = HEFF(I,J,bi,bj)            d_HSNWbyNEG(I,J)=MAX(-HSNOW(I,J,bi,bj),0._d 0)
282  # endif            HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)+d_HSNWbyNEG(I,J)
283              AREA(I,J,bi,bj)=MAX(AREA(I,J,bi,bj),0. _d 0)
284             ENDDO
285            ENDDO
286    
287    c 2) treat the case of very small area:
288    c
289    #ifdef ALLOW_AUTODIFF_TAMC
290    CADJ STORE area(:,:,bi,bj)  = comlev1_bibj,
291    CADJ &                        key = iicekey, byte = isbyte
292    #endif /* ALLOW_AUTODIFF_TAMC */
293            DO J=1,sNy
294             DO I=1,sNx
295              IF ((HEFF(i,j,bi,bj).GT.0).OR.(HSNOW(i,j,bi,bj).GT.0))
296         &     AREA(I,J,bi,bj)=MAX(AREA(I,J,bi,bj),areaMin)
297             ENDDO
298            ENDDO
299    
300    c 3) store regularized values of heff, hsnow, area at the onset of thermo.
301            DO J=1,sNy
302             DO I=1,sNx
303              HEFFpreTH(I,J)=HEFF(I,J,bi,bj)
304              HSNWpreTH(I,J)=HSNOW(I,J,bi,bj)
305              AREApreTH(I,J)=AREA(I,J,bi,bj)
306             ENDDO
307            ENDDO
308    
309    c 4) treat sea ice salinity pathological cases
310    #ifdef SEAICE_SALINITY
311    #ifdef ALLOW_AUTODIFF_TAMC
312    CADJ STORE hsalt(:,:,bi,bj)  = comlev1_bibj,
313    CADJ &                        key = iicekey, byte = isbyte
314    CADJ STORE heff(:,:,bi,bj)  = comlev1_bibj,
315    CADJ &                        key = iicekey, byte = isbyte
316    #endif /* ALLOW_AUTODIFF_TAMC */
317            DO J=1,sNy
318             DO I=1,sNx
319              IF ( (HSALT(I,J,bi,bj) .LT. 0.0).OR.
320         &         (HEFF(I,J,bi,bj) .EQ. 0.0)  ) THEN
321                 saltFluxAdjust(I,J) = - HEFFM(I,J,bi,bj) *
322         &            HSALT(I,J,bi,bj) / SEAICE_deltaTtherm
323                 HSALT(I,J,bi,bj) = 0.0 _d 0
324              ENDIF
325             ENDDO
326            ENDDO
327    #endif /* SEAICE_SALINITY */
328    
329    c 5) treat sea ice age pathological cases
330    c ...
331            
332    #else
333    #ifdef ALLOW_AUTODIFF_TAMC
334    CADJ STORE hsnow(:,:,bi,bj)  = comlev1_bibj,
335    CADJ &                        key = iicekey, byte = isbyte
336    #endif /* ALLOW_AUTODIFF_TAMC */
337            DO J=1,sNy
338             DO I=1,sNx
339              HEFFpreTH(I,J)=HEFFNM1(I,J,bi,bj)
340              HSNWpreTH(I,J)=HSNOW(I,J,bi,bj)
341              AREApreTH(I,J)=AREANM1(I,J,bi,bj)
342              d_HEFFbyNEG(I,J)=0. _d 0
343              d_HSNWbyNEG(I,J)=0. _d 0
344           ENDDO           ENDDO
345          ENDDO          ENDDO
346  #endif /* SEAICE_AGE */  #endif /* SEAICE_GROWTH_LEGACY */
347    
348    
349    c 4) COMPUTE ACTUAL ICE/SNOW THICKNESS; PUT MINIMUM ON HEFFACTUAL
350    c    TO REGULARIZE SEAICE_SOLVE4TEMP AND d_AREA COMPUTATIONS
351    c
352  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
353    CADJ STORE area(:,:,bi,bj)  = comlev1_bibj,
354    CADJ &                        key = iicekey, byte = isbyte
355  CADJ STORE heff(:,:,bi,bj)  = comlev1_bibj,  CADJ STORE heff(:,:,bi,bj)  = comlev1_bibj,
356  CADJ &                        key = iicekey, byte = isbyte  CADJ &                        key = iicekey, byte = isbyte
357  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,
# Line 281  CADJ &                        key = iice Line 359  CADJ &                        key = iice
359  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
360          DO J=1,sNy          DO J=1,sNy
361           DO I=1,sNx           DO I=1,sNx
362  C     COMPUTE ACTUAL ICE/SNOW THICKNESS AND PUT MINIMA            tmpscal1         = MAX(areaMin,AREApreTH(I,J))
363  C     TO REGULARIZE SEAICE_SOLVE4TEMP COMPUTATION            hsnowActual(I,J) = HSNWpreTH(I,J)/tmpscal1
364            areaSup          = MAX(areaMin,AREANm1(I,J,bi,bj))            tmpscal2         = HEFFpreTH(I,J)/tmpscal1
365            hsnowActual(I,J) = HSNOW(I,J,bi,bj)/areaSup            heffActual(I,J)  = MAX(tmpscal2,hiceMin)
           heffActual(I,J)  = HEFFNm1(I,J,bi,bj)/areaSup  
           heffActual(I,J)  = MAX(heffActual(I,J),hiceMin)  
366  C     Capping the actual ice thickness effectively enforces a  C     Capping the actual ice thickness effectively enforces a
367  C     minimum of heat flux through the ice and helps getting rid of  C     minimum of heat flux through the ice and helps getting rid of
368  C     very thick ice.  C     very thick ice.
# Line 297  cdm          heffActual(I,J)    = MIN(he Line 373  cdm          heffActual(I,J)    = MIN(he
373          ENDDO          ENDDO
374    
375    
376    
377  C determine available heat due to the atmosphere -- for open water  C determine available heat due to the atmosphere -- for open water
378  C ================================================================  C ================================================================
379    
# Line 400  C--  End loop over multi-categories Line 477  C--  End loop over multi-categories
477            DO J=1,sNy            DO J=1,sNy
478             DO I=1,sNx             DO I=1,sNx
479              DIAGarray(I,J) = maskC(I,J,kSurface,bi,bj) * (              DIAGarray(I,J) = maskC(I,J,kSurface,bi,bj) * (
480       &           a_QbyATM_cover(I,J) * areaNm1(I,J,bi,bj) +       &           a_QbyATM_cover(I,J) * AREApreTH(I,J) +
481       &           a_QbyATM_open(I,J) * ( ONE - areaNm1(I,J,bi,bj) ) )       &           a_QbyATM_open(I,J) * ( ONE - AREApreTH(I,J) ) )
482             ENDDO             ENDDO
483            ENDDO            ENDDO
484            CALL DIAGNOSTICS_FILL(DIAGarray,'SIatmQnt',0,1,3,bi,bj,myThid)            CALL DIAGNOSTICS_FILL(DIAGarray,'SIatmQnt',0,1,3,bi,bj,myThid)
# Line 413  c switch heat fluxes from W/m2 to 'effec Line 490  c switch heat fluxes from W/m2 to 'effec
490          DO J=1,sNy          DO J=1,sNy
491           DO I=1,sNx           DO I=1,sNx
492             a_QbyATM_cover(I,J) = a_QbyATM_cover(I,J)               a_QbyATM_cover(I,J) = a_QbyATM_cover(I,J)  
493       &         * convertQ2HI * areaNm1(I,J,bi,bj)       &         * convertQ2HI * AREApreTH(I,J)
494             a_QSWbyATM_cover(I,J) = a_QSWbyATM_cover(I,J)             a_QSWbyATM_cover(I,J) = a_QSWbyATM_cover(I,J)
495       &         * convertQ2HI * areaNm1(I,J,bi,bj)       &         * convertQ2HI * AREApreTH(I,J)
496             a_QbyATM_open(I,J) = a_QbyATM_open(I,J)             a_QbyATM_open(I,J) = a_QbyATM_open(I,J)
497       &         * convertQ2HI * ( ONE - areaNm1(I,J,bi,bj) )       &         * convertQ2HI * ( ONE - AREApreTH(I,J) )
498             a_QSWbyATM_open(I,J) = a_QSWbyATM_open(I,J)             a_QSWbyATM_open(I,J) = a_QSWbyATM_open(I,J)
499       &         * convertQ2HI * ( ONE - areaNm1(I,J,bi,bj) )       &         * convertQ2HI * ( ONE - AREApreTH(I,J) )
500  c and initialize r_QbyATM_cover  c and initialize r_QbyATM_cover
501             r_QbyATM_cover(I,J)=a_QbyATM_cover(I,J)             r_QbyATM_cover(I,J)=a_QbyATM_cover(I,J)
502           ENDDO           ENDDO
# Line 432  C ====================================== Line 509  C ======================================
509          DO J=1,sNy          DO J=1,sNy
510           DO I=1,sNx           DO I=1,sNx
511            IF (a_QbyATM_cover(I,J).LT.ZERO.AND.            IF (a_QbyATM_cover(I,J).LT.ZERO.AND.
512       &        AREANm1(I,J,bi,bj).GT.ZERO) THEN       &        AREApreTH(I,J).GT.ZERO) THEN
513              FRWfromSNW(I,J)=2              FRWfromSNW(I,J)=2
514            ELSE            ELSE
515              FRWfromSNW(I,J)=1              FRWfromSNW(I,J)=1
# Line 540  C ====================================== Line 617  C ======================================
617  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
618  CADJ STORE heff(:,:,bi,bj)   = comlev1_bibj,  CADJ STORE heff(:,:,bi,bj)   = comlev1_bibj,
619  CADJ &                         key = iicekey, byte = isbyte  CADJ &                         key = iicekey, byte = isbyte
620    CADJ STORE area(:,:,bi,bj)   = comlev1_bibj,
621    CADJ &                         key = iicekey, byte = isbyte
622  CADJ STORE r_QbyATM_cover(:,:)         = comlev1_bibj,  CADJ STORE r_QbyATM_cover(:,:)         = comlev1_bibj,
623  CADJ &                         key = iicekey, byte = isbyte  CADJ &                         key = iicekey, byte = isbyte
624  CADJ STORE a_QbyATM_cover(:,:)        = comlev1_bibj,  CADJ STORE a_QbyATM_cover(:,:)        = comlev1_bibj,
# Line 567  C gain of new ice over open water (>0 by Line 646  C gain of new ice over open water (>0 by
646  c compute cover fraction tendency  c compute cover fraction tendency
647            IF ( YC(I,J,bi,bj) .LT. ZERO ) THEN            IF ( YC(I,J,bi,bj) .LT. ZERO ) THEN
648             d_AREAbyATM(I,J)=tmpscal4/HO_south             d_AREAbyATM(I,J)=tmpscal4/HO_south
649       &          +HALF*tmpscal3*AREANm1(I,J,bi,bj)  #ifndef SEAICE_GROWTH_LEGACY
650         &          +HALF*tmpscal3/heffActual(I,J)
651    #else
652         &          +HALF*tmpscal3*AREApreTH(I,J)
653       &          /(HEFF(I,J,bi,bj)+.00001 _d 0)       &          /(HEFF(I,J,bi,bj)+.00001 _d 0)
654    #endif
655            ELSE            ELSE
656             d_AREAbyATM(I,J)=tmpscal4/HO             d_AREAbyATM(I,J)=tmpscal4/HO
657       &          +HALF*tmpscal3*AREANm1(I,J,bi,bj)  #ifndef SEAICE_GROWTH_LEGACY
658         &          +HALF*tmpscal3/heffActual(I,J)
659    #else
660         &          +HALF*tmpscal3*AREApreTH(I,J)
661       &          /(HEFF(I,J,bi,bj)+.00001 _d 0)       &          /(HEFF(I,J,bi,bj)+.00001 _d 0)
662    #endif
663            ENDIF            ENDIF
664  c apply tendency  c apply tendency
665             AREA(I,J,bi,bj)=AREA(I,J,bi,bj)+d_AREAbyATM(I,J)             AREA(I,J,bi,bj)=max(0. _d 0 , min( 1. _d 0,
666         &                     AREA(I,J,bi,bj)+d_AREAbyATM(I,J) ) )
667           ENDDO           ENDDO
668          ENDDO          ENDDO
669    
# Line 626  C ====================================== Line 714  C ======================================
714  C           add precip as snow  C           add precip as snow
715              d_HFRWbyRAIN(I,J)=0. _d 0              d_HFRWbyRAIN(I,J)=0. _d 0
716              d_HSNWbyRAIN(I,J)=convertPRECIP2HI*ICE2SNOW*              d_HSNWbyRAIN(I,J)=convertPRECIP2HI*ICE2SNOW*
717       &            PRECIP(I,J,bi,bj)*AREANm1(I,J,bi,bj)       &            PRECIP(I,J,bi,bj)*AREApreTH(I,J)
718            ELSE            ELSE
719  c           add precip to the fresh water bucket  c           add precip to the fresh water bucket
720              d_HFRWbyRAIN(I,J)=-convertPRECIP2HI*              d_HFRWbyRAIN(I,J)=-convertPRECIP2HI*
721       &            PRECIP(I,J,bi,bj)*AREANm1(I,J,bi,bj)       &            PRECIP(I,J,bi,bj)*AREApreTH(I,J)
722              d_HSNWbyRAIN(I,J)=0. _d 0              d_HSNWbyRAIN(I,J)=0. _d 0
723            ENDIF            ENDIF
724  c apply tendency  c apply tendency
# Line 668  cgf heat and water conservation: ok Line 756  cgf heat and water conservation: ok
756  #endif /* SEAICE_EXCLUDE_FOR_EXACT_AD_TESTING */  #endif /* SEAICE_EXCLUDE_FOR_EXACT_AD_TESTING */
757  cph)  cph)
758    
759    #ifndef SEAICE_GROWTH_LEGACY
760    C convert snow to ice if submerged.
761    C =================================
762    
763    #ifdef ALLOW_SEAICE_FLOODING
764            IF ( SEAICEuseFlooding ) THEN
765             DO J=1,sNy
766              DO I=1,sNx
767               hDraft     = (HSNOW(I,J,bi,bj)*SEAICE_rhoSnow
768         &              +HEFF(I,J,bi,bj)*SEAICE_rhoIce)/rhoConst
769    cgf - use of 1000 instead of rhoConst (or rhoConstFresh?)
770               tmpscal1 = MAX( 0. _d 0, hDraft - HEFF(I,J,bi,bj))
771               d_HEFFbyFLOODING(I,J)=tmpscal1
772    c apply tendency
773               HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj)+d_HEFFbyFLOODING(I,J)
774               HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)-
775         &                           d_HEFFbyFLOODING(I,J)*ICE2SNOW  
776    cgf heat and water conservation: ok -- since ice and snow
777    cgf are treated as having a consistent latent heat of fusion
778              ENDDO
779             ENDDO
780    #ifdef ALLOW_DIAGNOSTICS
781             IF ( useDiagnostics ) THEN
782              IF ( DIAGNOSTICS_IS_ON('SIsnwice',myThid) ) THEN
783    C     turn tmparr1 into a rate
784               DO J=1,sNy
785                DO I=1,sNx
786                 tmparr1(I,J) = d_HEFFbyFLOODING(I,J)/SEAICE_deltaTtherm
787                ENDDO
788               ENDDO
789               CALL DIAGNOSTICS_FILL(tmparr1,'SIsnwice',0,1,3,bi,bj,myThid)
790              ENDIF
791             ENDIF
792    #endif /* ALLOW_DIAGNOSTICS */
793            ENDIF
794    #endif /* ALLOW_SEAICE_FLOODING */
795    #endif /* SEAICE_GROWTH_LEGACY */
796    
797    
798    
799  C compute net fresh water flux leaving/entering  C compute net fresh water flux leaving/entering
800  C the ocean, accounting for fresh/salt water stocks.  C the ocean, accounting for fresh/salt water stocks.
# Line 690  CADJ &                        key = iice Line 817  CADJ &                        key = iice
817       &            +d_HSNWbyOCNonSNW(I,J)/ICE2SNOW       &            +d_HSNWbyOCNonSNW(I,J)/ICE2SNOW
818       &            +d_HEFFbyOCNonICE(I,J)       &            +d_HEFFbyOCNonICE(I,J)
819       &            +d_HEFFbyATMonOCN(I,J)       &            +d_HEFFbyATMonOCN(I,J)
820         &            +d_HEFFbyNEG(I,J)
821         &            +d_HSNWbyNEG(I,J)/ICE2SNOW
822  C NOW GET FRESH WATER FLUX  C NOW GET FRESH WATER FLUX
823            EmPmR(I,J,bi,bj)  = maskC(I,J,kSurface,bi,bj)*(            EmPmR(I,J,bi,bj)  = maskC(I,J,kSurface,bi,bj)*(
824       &         ( EVAP(I,J,bi,bj)-PRECIP(I,J,bi,bj) )       &         ( EVAP(I,J,bi,bj)-PRECIP(I,J,bi,bj) )
825       &         * ( ONE - AREANm1(I,J,bi,bj) )       &         * ( ONE - AREApreTH(I,J) )
826  #ifdef ALLOW_RUNOFF  #ifdef ALLOW_RUNOFF
827       &         - RUNOFF(I,J,bi,bj)       &         - RUNOFF(I,J,bi,bj)
828  #endif  #endif
# Line 711  cgf heat and water conservation: ok -- s Line 840  cgf heat and water conservation: ok -- s
840              DIAGarray(I,J) = maskC(I,J,kSurface,bi,bj)*(              DIAGarray(I,J) = maskC(I,J,kSurface,bi,bj)*(
841       &           PRECIP(I,J,bi,bj)       &           PRECIP(I,J,bi,bj)
842       &           - EVAP(I,J,bi,bj)       &           - EVAP(I,J,bi,bj)
843       &           *( ONE - AREANm1(I,J,bi,bj) )       &           *( ONE - AREApreTH(I,J) )
844       &           + RUNOFF(I,J,bi,bj)       &           + RUNOFF(I,J,bi,bj)
845       &           )*rhoConstFresh       &           )*rhoConstFresh
846             ENDDO             ENDDO
# Line 726  cgf heat and water conservation: ok -- s Line 855  cgf heat and water conservation: ok -- s
855            frWtrAtm(I,J,bi,bj) = maskC(I,J,kSurface,bi,bj)*(            frWtrAtm(I,J,bi,bj) = maskC(I,J,kSurface,bi,bj)*(
856       &         PRECIP(I,J,bi,bj)       &         PRECIP(I,J,bi,bj)
857       &         - EVAP(I,J,bi,bj)       &         - EVAP(I,J,bi,bj)
858       &         *( ONE - AREANm1(I,J,bi,bj) )       &         *( ONE - AREApreTH(I,J) )
859       &         + RUNOFF(I,J,bi,bj)       &         + RUNOFF(I,J,bi,bj)
860       &         )*rhoConstFresh       &         )*rhoConstFresh
861           ENDDO           ENDDO
# Line 739  C ====================================== Line 868  C ======================================
868    
869  #ifdef SEAICE_SALINITY  #ifdef SEAICE_SALINITY
870    
871    #ifdef SEAICE_GROWTH_LEGACY
872          DO J=1,sNy          DO J=1,sNy
873           DO I=1,sNx           DO I=1,sNx
874  C set HSALT = 0 if HSALT < 0 and compute salt to remove from ocean  C set HSALT = 0 if HSALT < 0 and compute salt to remove from ocean
# Line 749  C set HSALT = 0 if HSALT < 0 and compute Line 879  C set HSALT = 0 if HSALT < 0 and compute
879            ENDIF            ENDIF
880           ENDDO           ENDDO
881          ENDDO          ENDDO
882    #endif /* SEAICE_GROWTH_LEGACY */
883    
884  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
885  CADJ STORE hsalt(:,:,bi,bj) = comlev1_bibj,  CADJ STORE hsalt(:,:,bi,bj) = comlev1_bibj,
# Line 757  CADJ &                        key = iice Line 888  CADJ &                        key = iice
888    
889          DO J=1,sNy          DO J=1,sNy
890           DO I=1,sNx           DO I=1,sNx
891    c sum up the terms that affect the salt content of the ice pack
892           tmpscal1=d_HEFFbyOCNonICE(I,J)+d_HEFFbyATMonOCN(I,J)           tmpscal1=d_HEFFbyOCNonICE(I,J)+d_HEFFbyATMonOCN(I,J)
893    c recompute HEFF before thermodyncamic updates (which is not AREApreTH in legacy code)
894             tmpscal2=HEFF(I,J,bi,bj)-tmpscal1-d_HEFFbyFLOODING(I,J)
895  C tmpscal1 > 0 : m of sea ice that is created  C tmpscal1 > 0 : m of sea ice that is created
896            IF ( tmpscal1 .GE. 0.0 ) THEN            IF ( tmpscal1 .GE. 0.0 ) THEN
897               saltFlux(I,J,bi,bj) =               saltFlux(I,J,bi,bj) =
898       &            HEFFM(I,J,bi,bj)*tmpscal1*       &            HEFFM(I,J,bi,bj)/SEAICE_deltaTtherm
899       &            ICE2WATR*rhoConstFresh*SEAICE_salinity*       &            *SEAICE_salinity*salt(I,j,kSurface,bi,bj)
900       &            salt(I,j,kSurface,bi,bj)/SEAICE_deltaTtherm       &            *tmpscal1*ICE2WATR*rhoConstFresh
901  #ifdef ALLOW_SALT_PLUME  #ifdef ALLOW_SALT_PLUME
902  C saltPlumeFlux is defined only during freezing:  C saltPlumeFlux is defined only during freezing:
903               saltPlumeFlux(I,J,bi,bj)=               saltPlumeFlux(I,J,bi,bj)=
904       &            HEFFM(I,J,bi,bj)*tmpscal1*       &            HEFFM(I,J,bi,bj)/SEAICE_deltaTtherm
905       &            ICE2WATR*rhoConstFresh*(1-SEAICE_salinity)*       &            *(1-SEAICE_salinity)*salt(I,j,kSurface,bi,bj)
906       &            salt(I,j,kSurface,bi,bj)/SEAICE_deltaTtherm       &            *tmpscal1*ICE2WATR*rhoConstFresh
907  C if SaltPlumeSouthernOcean=.FALSE. turn off salt plume in Southern Ocean  C if SaltPlumeSouthernOcean=.FALSE. turn off salt plume in Southern Ocean
908               IF ( .NOT. SaltPlumeSouthernOcean ) THEN               IF ( .NOT. SaltPlumeSouthernOcean ) THEN
909                IF ( YC(I,J,bi,bj) .LT. 0.0 _d 0 )                IF ( YC(I,J,bi,bj) .LT. 0.0 _d 0 )
# Line 780  C if SaltPlumeSouthernOcean=.FALSE. turn Line 914  C if SaltPlumeSouthernOcean=.FALSE. turn
914  C tmpscal1 < 0 : m of sea ice that is melted  C tmpscal1 < 0 : m of sea ice that is melted
915            ELSE            ELSE
916               saltFlux(I,J,bi,bj) =               saltFlux(I,J,bi,bj) =
917       &            HEFFM(I,J,bi,bj)*tmpscal1*       &         HEFFM(I,J,bi,bj)/SEAICE_deltaTtherm
918       &            HSALT(I,J,bi,bj)/       &         *HSALT(I,J,bi,bj)
919       &            (HEFF(I,J,bi,bj)-tmpscal1)/       &         *tmpscal1/tmpscal2
      &            SEAICE_deltaTtherm  
920  #ifdef ALLOW_SALT_PLUME  #ifdef ALLOW_SALT_PLUME
921               saltPlumeFlux(i,j,bi,bj) = 0.0 _d 0               saltPlumeFlux(i,j,bi,bj) = 0.0 _d 0
922  #endif /* ALLOW_SALT_PLUME */  #endif /* ALLOW_SALT_PLUME */
# Line 793  C update HSALT based on surface saltFlux Line 926  C update HSALT based on surface saltFlux
926       &         saltFlux(I,J,bi,bj) * SEAICE_deltaTtherm       &         saltFlux(I,J,bi,bj) * SEAICE_deltaTtherm
927            saltFlux(I,J,bi,bj) =            saltFlux(I,J,bi,bj) =
928       &         saltFlux(I,J,bi,bj) + saltFluxAdjust(I,J)       &         saltFlux(I,J,bi,bj) + saltFluxAdjust(I,J)
929    #ifdef SEAICE_GROWTH_LEGACY
930  C set HSALT = 0 if HEFF = 0 and compute salt to dump into ocean  C set HSALT = 0 if HEFF = 0 and compute salt to dump into ocean
931            IF ( HEFF(I,J,bi,bj) .EQ. 0.0 ) THEN            IF ( HEFF(I,J,bi,bj) .EQ. 0.0 ) THEN
932               saltFlux(I,J,bi,bj) = saltFlux(I,J,bi,bj) -               saltFlux(I,J,bi,bj) = saltFlux(I,J,bi,bj) -
# Line 803  C set HSALT = 0 if HEFF = 0 and compute Line 937  C set HSALT = 0 if HEFF = 0 and compute
937               saltPlumeFlux(i,j,bi,bj) = 0.0 _d 0               saltPlumeFlux(i,j,bi,bj) = 0.0 _d 0
938  #endif /* ALLOW_SALT_PLUME */  #endif /* ALLOW_SALT_PLUME */
939            ENDIF            ENDIF
940    #endif /* SEAICE_GROWTH_LEGACY */
941           ENDDO           ENDDO
942          ENDDO          ENDDO
943  #endif /* SEAICE_SALINITY */  #endif /* SEAICE_SALINITY */
# Line 818  C ====================================== Line 953  C ======================================
953            QNET(I,J,bi,bj) = r_QbyATM_cover(I,J) + a_QbyATM_open(I,J)            QNET(I,J,bi,bj) = r_QbyATM_cover(I,J) + a_QbyATM_open(I,J)
954  C account for contribution due to ocean-ice interaction.  C account for contribution due to ocean-ice interaction.
955       &         - ( d_HEFFbyOCNonICE(I,J) +       &         - ( d_HEFFbyOCNonICE(I,J) +
956       &             d_HSNWbyOCNonSNW(I,J)/ICE2SNOW )       &             d_HSNWbyOCNonSNW(I,J)/ICE2SNOW +
957         &             d_HEFFbyNEG(I,J) +
958         &             d_HSNWbyNEG(I,J)/ICE2SNOW )
959       &         * maskC(I,J,kSurface,bi,bj)       &         * maskC(I,J,kSurface,bi,bj)
960            QSW(I,J,bi,bj)  = a_QSWbyATM_cover(I,J) + a_QSWbyATM_open(I,J)            QSW(I,J,bi,bj)  = a_QSWbyATM_cover(I,J) + a_QSWbyATM_open(I,J)
961           ENDDO           ENDDO
# Line 829  C account for contribution due to ocean- Line 966  C account for contribution due to ocean-
966           IF ( DIAGNOSTICS_IS_ON('SIqneto ',myThid) ) THEN           IF ( DIAGNOSTICS_IS_ON('SIqneto ',myThid) ) THEN
967            DO J=1,sNy            DO J=1,sNy
968             DO I=1,sNx             DO I=1,sNx
969              DIAGarray(I,J) = a_QbyATM_open(I,J) *              DIAGarray(I,J) = a_QbyATM_open(I,J) * (ONE-AREApreTH(I,J))
      &         (ONE-areaNm1(I,J,bi,bj))  
970             ENDDO             ENDDO
971            ENDDO            ENDDO
972            CALL DIAGNOSTICS_FILL(DIAGarray,'SIqneto ',0,1,3,bi,bj,myThid)            CALL DIAGNOSTICS_FILL(DIAGarray,'SIqneto ',0,1,3,bi,bj,myThid)
# Line 838  C account for contribution due to ocean- Line 974  C account for contribution due to ocean-
974           IF ( DIAGNOSTICS_IS_ON('SIqneti ',myThid) ) THEN           IF ( DIAGNOSTICS_IS_ON('SIqneti ',myThid) ) THEN
975            DO J=1,sNy            DO J=1,sNy
976             DO I=1,sNx             DO I=1,sNx
977              DIAGarray(I,J) = r_QbyATM_cover(I,J) * areaNm1(I,J,bi,bj)              DIAGarray(I,J) = r_QbyATM_cover(I,J) * AREApreTH(I,J)
978             ENDDO             ENDDO
979            ENDDO            ENDDO
980            CALL DIAGNOSTICS_FILL(DIAGarray,'SIqneti ',0,1,3,bi,bj,myThid)            CALL DIAGNOSTICS_FILL(DIAGarray,'SIqneti ',0,1,3,bi,bj,myThid)
# Line 847  C account for contribution due to ocean- Line 983  C account for contribution due to ocean-
983  #endif  #endif
984    
985    
986    #ifdef SEAICE_GROWTH_LEGACY
987    
988  C treat values of ice cover fraction oustide  C treat values of ice cover fraction oustide
989  C the [0 1] range, and other such issues.  C the [0 1] range, and other such issues.
990  C ===========================================  C ===========================================
# Line 903  CADJ &                         key = iic Line 1041  CADJ &                         key = iic
1041  C     use (abuse) tmparr1 to diagnose the total thermodynamic growth rate  C     use (abuse) tmparr1 to diagnose the total thermodynamic growth rate
1042            DO J=1,sNy            DO J=1,sNy
1043             DO I=1,sNx             DO I=1,sNx
1044              tmparr1(I,J) = (HEFF(I,J,bi,bj)-HEFFNm1(I,J,bi,bj))              tmparr1(I,J) = (HEFF(I,J,bi,bj)-HEFFpreTH(I,J))
1045       &           /SEAICE_deltaTtherm       &           /SEAICE_deltaTtherm
1046             ENDDO             ENDDO
1047            ENDDO            ENDDO
# Line 949  C     turn tmparr1 into a rate Line 1087  C     turn tmparr1 into a rate
1087          ENDIF          ENDIF
1088  #endif /* ALLOW_SEAICE_FLOODING */  #endif /* ALLOW_SEAICE_FLOODING */
1089    
1090    #endif /* SEAICE_GROWTH_LEGACY */
1091    
1092  C Sea Ice Load on the sea surface.  C Sea Ice Load on the sea surface.
1093  C =================================  C =================================
# Line 974  C                 b) melting: ice fracti Line 1113  C                 b) melting: ice fracti
1113          DO J=1,sNy          DO J=1,sNy
1114           DO I=1,sNx           DO I=1,sNx
1115            IF ( AREA(I,J,bi,bj) .GT. 0.15 ) THEN            IF ( AREA(I,J,bi,bj) .GT. 0.15 ) THEN
1116             IF ( AREA(i,j,bi,bj) .LT. old_AREA(i,j) ) THEN             IF ( AREA(i,j,bi,bj) .LT. AREApreTH(i,j) ) THEN
1117  C--   scale effective ice-age to account for ice-age sink associated with melting  C--   scale effective ice-age to account for ice-age sink associated with melting
1118              IceAge(i,j,bi,bj) = IceAge(i,j,bi,bj)              IceAge(i,j,bi,bj) = IceAge(i,j,bi,bj)
1119       &         *AREA(i,j,bi,bj)/old_AREA(i,j)       &         *AREA(i,j,bi,bj)/AREApreTH(i,j)
1120             ENDIF             ENDIF
1121  C--   account for aging:  C--   account for aging:
1122             IceAge(i,j,bi,bj) = IceAge(i,j,bi,bj)             IceAge(i,j,bi,bj) = IceAge(i,j,bi,bj)
# Line 994  C                 b) melting: ice volume Line 1133  C                 b) melting: ice volume
1133          DO J=1,sNy          DO J=1,sNy
1134           DO I=1,sNx           DO I=1,sNx
1135  C--   compute actual age from effective age:  C--   compute actual age from effective age:
1136            IF (OLD_AREA(i,j).GT.0. _d 0) THEN            IF (AREApreTH(i,j).GT.0. _d 0) THEN
1137             age_actual=IceAge(i,j,bi,bj)/OLD_AREA(i,j)             tmpscal1=IceAge(i,j,bi,bj)/AREApreTH(i,j)
1138            ELSE            ELSE
1139             age_actual=0. _d 0             tmpscal1=0. _d 0
1140            ENDIF                    ENDIF        
1141            IF ( (OLD_HEFF(i,j).LT.HEFF(i,j,bi,bj)).AND.            IF ( (HEFFpreTH(i,j).LT.HEFF(i,j,bi,bj)).AND.
1142       &         (AREA(i,j,bi,bj).GT.0.15) ) THEN       &         (AREA(i,j,bi,bj).GT.0.15) ) THEN
1143             age_actual=age_actual*OLD_HEFF(i,j)/             tmpscal2=tmpscal1*HEFFpreTH(i,j)/
1144       &          HEFF(i,j,bi,bj)+SEAICE_deltaTtherm       &          HEFF(i,j,bi,bj)+SEAICE_deltaTtherm
1145            ELSEIF (AREA(i,j,bi,bj).LE.0.15) THEN            ELSEIF (AREA(i,j,bi,bj).LE.0.15) THEN
1146             age_actual=0. _d 0             tmpscal2=0. _d 0
1147            ELSE            ELSE
1148             age_actual=age_actual+SEAICE_deltaTtherm             tmpscal2=tmpscal1+SEAICE_deltaTtherm
1149            ENDIF            ENDIF
1150  C--   re-scale to effective age:  C--   re-scale to effective age:
1151            IceAge(i,j,bi,bj) = age_actual*AREA(i,j,bi,bj)            IceAge(i,j,bi,bj) = tmpscal2*AREA(i,j,bi,bj)
1152           ENDDO           ENDDO
1153          ENDDO          ENDDO
1154  # endif /* SEAICE_AGE_VOL */  # endif /* SEAICE_AGE_VOL */
Legend:
Removed from v.1.86  
changed lines
  Added in v.1.87
  ViewVC Help
Powered by ViewVC 1.1.22