/[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.60 by mlosch, Mon Jun 22 15:57:45 2009 UTC revision 1.61 by mlosch, Wed Jun 24 08:01:43 2009 UTC
# Line 147  C Line 147  C
147  C     initialise a few fields  C     initialise a few fields
148  C  C
149  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
150  CADJ STORE area(:,:,:,bi,bj) = comlev1_bibj,  CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
151  CADJ &                         key = iicekey, byte = isbyte  CADJ &                       key = iicekey, byte = isbyte
152  CADJ STORE qnet(:,:,bi,bj)   = comlev1_bibj,  CADJ STORE qnet(:,:,bi,bj) = comlev1_bibj,
153  CADJ &                         key = iicekey, byte = isbyte  CADJ &                       key = iicekey, byte = isbyte
154  CADJ STORE qsw(:,:,bi,bj)    = comlev1_bibj,  CADJ STORE qsw(:,:,bi,bj)  = comlev1_bibj,
155  CADJ &                         key = iicekey, byte = isbyte  CADJ &                       key = iicekey, byte = isbyte
156  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
157          DO J=1,sNy          DO J=1,sNy
158           DO I=1,sNx           DO I=1,sNx
# Line 180  CADJ &                         key = iic Line 180  CADJ &                         key = iic
180           ENDDO           ENDDO
181          ENDDO          ENDDO
182  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
183  CADJ STORE heff(:,:,:,bi,bj) = comlev1_bibj,  CADJ STORE heff(:,:,bi,bj)  = comlev1_bibj,
184  CADJ &                         key = iicekey, byte = isbyte  CADJ &                        key = iicekey, byte = isbyte
185  CADJ STORE hsnow(:,:,bi,bj)  = comlev1_bibj,  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,
186  CADJ &                         key = iicekey, byte = isbyte  CADJ &                        key = iicekey, byte = isbyte
187  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
188          DO J=1,sNy          DO J=1,sNy
189           DO I=1,sNx           DO I=1,sNx
# Line 259  C     Compute relative wind speed over s Line 259  C     Compute relative wind speed over s
259            DO I=1,sNx            DO I=1,sNx
260             SPEED_SQ =             SPEED_SQ =
261       &          (uWind(I,J,bi,bj)       &          (uWind(I,J,bi,bj)
262       &          +0.5 _d 0*(uVel(i,j,1,bi,bj)+uVel(i+1,j,1,bi,bj))       &          +0.5 _d 0*(uVel(i,j,kSurface,bi,bj)
263       &          -0.5 _d 0*(uice(i,j,1,bi,bj)+uice(i+1,j,1,bi,bj)))**2       &                    +uVel(i+1,j,kSurface,bi,bj))
264         &          -0.5 _d 0*(uice(i,j,bi,bj)+uice(i+1,j,bi,bj)))**2
265       &          +(vWind(I,J,bi,bj)       &          +(vWind(I,J,bi,bj)
266       &          +0.5 _d 0*(vVel(i,j,1,bi,bj)+vVel(i,j+1,1,bi,bj))       &          +0.5 _d 0*(vVel(i,j,kSurface,bi,bj)
267       &          -0.5 _d 0*(vice(i,j,1,bi,bj)+vice(i,j+1,1,bi,bj)))**2       &                    +vVel(i,j+1,kSurface,bi,bj))
268         &          -0.5 _d 0*(vice(i,j,bi,bj)+vice(i,j+1,bi,bj)))**2
269             IF ( SPEED_SQ .LE. SEAICE_EPS_SQ ) THEN             IF ( SPEED_SQ .LE. SEAICE_EPS_SQ ) THEN
270               UG(I,J)=SEAICE_EPS               UG(I,J)=SEAICE_EPS
271             ELSE             ELSE
# Line 316  C--  End loop over multi-categories Line 318  C--  End loop over multi-categories
318            DO J=1,sNy            DO J=1,sNy
319             DO I=1,sNx             DO I=1,sNx
320              DIAGarray(I,J) = maskC(I,J,kSurface,bi,bj) * (              DIAGarray(I,J) = maskC(I,J,kSurface,bi,bj) * (
321       &           FICE(I,J) * AREA(I,J,2,bi,bj) +       &           FICE(I,J) * areaNm1(I,J,bi,bj) +
322       &           QNETO(I,J) * ( ONE - AREA(I,J,2,bi,bj) ) )       &           QNETO(I,J) * ( ONE - areaNm1(I,J,bi,bj) ) )
323             ENDDO             ENDDO
324            ENDDO            ENDDO
325            CALL DIAGNOSTICS_FILL(DIAGarray,'SIatmQnt',0,1,3,bi,bj,myThid)            CALL DIAGNOSTICS_FILL(DIAGarray,'SIatmQnt',0,1,3,bi,bj,myThid)
# Line 326  C--  End loop over multi-categories Line 328  C--  End loop over multi-categories
328  #endif  #endif
329    
330  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
331  CADJ STORE theta(:,:,:,bi,bj)= comlev1_bibj,  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj,
332  CADJ &                         key = iicekey, byte = isbyte  CADJ &                          key = iicekey, byte = isbyte
333  CADJ STORE heff(:,:,:,bi,bj) = comlev1_bibj,  CADJ STORE heff(:,:,bi,bj)    = comlev1_bibj,
334  CADJ &                         key = iicekey, byte = isbyte  CADJ &                          key = iicekey, byte = isbyte
335  #endif  #endif
336  C  C
337  C--   compute and apply ice growth due to oceanic heat flux from below  C--   compute and apply ice growth due to oceanic heat flux from below
# Line 364  C     Positive YNEG/availHeat leads to i Line 366  C     Positive YNEG/availHeat leads to i
366             availHeat(i,j) = 0.             availHeat(i,j) = 0.
367            ENDIF            ENDIF
368  C     local copy of old effective ice thickness  C     local copy of old effective ice thickness
369            hEffOld(I,J)      = HEFF(I,J,1,bi,bj)            hEffOld(I,J)      = HEFF(I,J,bi,bj)
370  C     Melt (availHeat>0) or create (availHeat<0) sea ice  C     Melt (availHeat>0) or create (availHeat<0) sea ice
371            HEFF(I,J,1,bi,bj) = MAX(ZERO,HEFF(I,J,1,bi,bj)-availHeat(I,J))            HEFF(I,J,bi,bj) = MAX(ZERO,HEFF(I,J,bi,bj)-availHeat(I,J))
372  C      C    
373            YNEG(I,J,bi,bj)   = hEffOld(I,J)    - HEFF(I,J,1,bi,bj)            YNEG(I,J,bi,bj)   = hEffOld(I,J)    - HEFF(I,J,bi,bj)
374  C      C    
375            saltWtrIce(I,J,bi,bj)   = saltWtrIce(I,J,bi,bj)            saltWtrIce(I,J,bi,bj)   = saltWtrIce(I,J,bi,bj)
376       &                              - YNEG(I,J,bi,bj)       &                              - YNEG(I,J,bi,bj)
# Line 395  cphCADJ STORE hsnow  = comlev1, key = ik Line 397  cphCADJ STORE hsnow  = comlev1, key = ik
397  cph)  cph)
398  c  c
399  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
400  CADJ STORE area(:,:,:,bi,bj) = comlev1_bibj,  CADJ STORE area(:,:,bi,bj)   = comlev1_bibj,
401  CADJ &                         key = iicekey, byte = isbyte  CADJ &                         key = iicekey, byte = isbyte
402  CADJ STORE hsnow(:,:,bi,bj)  = comlev1_bibj,  CADJ STORE hsnow(:,:,bi,bj)  = comlev1_bibj,
403  CADJ &                         key = iicekey, byte = isbyte  CADJ &                         key = iicekey, byte = isbyte
# Line 460  C     now get cell averaged growth rate Line 462  C     now get cell averaged growth rate
462          ENDDO          ENDDO
463  cph(  cph(
464  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
465  CADJ STORE heff(:,:,:,bi,bj) = comlev1_bibj,  CADJ STORE heff(:,:,bi,bj)   = comlev1_bibj,
466  CADJ &                         key = iicekey, byte = isbyte  CADJ &                         key = iicekey, byte = isbyte
467  CADJ STORE fice(:,:)         = comlev1_bibj,  CADJ STORE fice(:,:)         = comlev1_bibj,
468  CADJ &                         key = iicekey, byte = isbyte  CADJ &                         key = iicekey, byte = isbyte
# Line 483  C     negative growth in meters of ice ( Line 485  C     negative growth in meters of ice (
485            growthNeg  = -SEAICE_deltaTtherm*FHEFF(I,J)*recip_QI            growthNeg  = -SEAICE_deltaTtherm*FHEFF(I,J)*recip_QI
486  C     negative growth must not exceed effective ice thickness (=volume)  C     negative growth must not exceed effective ice thickness (=volume)
487  C     (that is, cannot melt more than all the ice)  C     (that is, cannot melt more than all the ice)
488            growthHEFF = -ONE*MIN(HEFF(I,J,1,bi,bj),growthNeg)            growthHEFF = -ONE*MIN(HEFF(I,J,bi,bj),growthNeg)
489  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
490            DIAGarray(I,J) = growthHEFF            DIAGarray(I,J) = growthHEFF
491  #endif  #endif
# Line 495  CML   removed these loops and moved TAMC Line 497  CML   removed these loops and moved TAMC
497  CML         ENDDO  CML         ENDDO
498  CML        ENDDO  CML        ENDDO
499  CML#ifdef ALLOW_AUTODIFF_TAMC  CML#ifdef ALLOW_AUTODIFF_TAMC
500  CMLCADJ STORE area(:,:,:,bi,bj) = comlev1_bibj,  CMLCADJ STORE area(:,:,bi,bj) = comlev1_bibj,
501  CMLCADJ &                         key = iicekey, byte = isbyte  CMLCADJ &                       key = iicekey, byte = isbyte
502  CML#endif  CML#endif
503  CML        DO J=1,sNy  CML        DO J=1,sNy
504  CML         DO I=1,sNx  CML         DO I=1,sNx
505  C     Here we finally compute the new AREA  C     Here we finally compute the new AREA
506            IF ( YC(I,J,bi,bj) .LT. ZERO ) THEN            IF ( YC(I,J,bi,bj) .LT. ZERO ) THEN
507             AREA(I,J,1,bi,bj)=AREA(I,J,1,bi,bj)+             AREA(I,J,bi,bj)=AREA(I,J,bi,bj)+
508       &          (ONE-AREANm1(I,J,bi,bj))*GAREA(I,J)/HO_south       &          (ONE-AREANm1(I,J,bi,bj))*GAREA(I,J)/HO_south
509       &          +HALF*HCORR(I,J)*AREANm1(I,J,bi,bj)       &          +HALF*HCORR(I,J)*AREANm1(I,J,bi,bj)
510       &          /(HEFF(I,J,1,bi,bj)+.00001 _d 0)       &          /(HEFF(I,J,bi,bj)+.00001 _d 0)
511            ELSE            ELSE
512             AREA(I,J,1,bi,bj)=AREA(I,J,1,bi,bj)+             AREA(I,J,bi,bj)=AREA(I,J,bi,bj)+
513       &          (ONE-AREANm1(I,J,bi,bj))*GAREA(I,J)/HO       &          (ONE-AREANm1(I,J,bi,bj))*GAREA(I,J)/HO
514       &          +HALF*HCORR(I,J)*AREANm1(I,J,bi,bj)       &          +HALF*HCORR(I,J)*AREANm1(I,J,bi,bj)
515       &          /(HEFF(I,J,1,bi,bj)+.00001 _d 0)       &          /(HEFF(I,J,bi,bj)+.00001 _d 0)
516            ENDIF            ENDIF
517           ENDDO           ENDDO
518          ENDDO          ENDDO
# Line 524  C     Here we finally compute the new AR Line 526  C     Here we finally compute the new AR
526  #endif  #endif
527    
528  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
529  CADJ STORE area(:,:,:,bi,bj) = comlev1_bibj,  CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
530  CADJ &                         key = iicekey, byte = isbyte  CADJ &                       key = iicekey, byte = isbyte
531  #endif  #endif
532  C      C    
533  C     now update (freeze or melt) HEFF  C     now update (freeze or melt) HEFF
# Line 533  C Line 535  C
535          DO J=1,sNy          DO J=1,sNy
536           DO I=1,sNx           DO I=1,sNx
537  C     negative growth (>0 for melting) of existing ice in meters  C     negative growth (>0 for melting) of existing ice in meters
538            growthNeg        = -SEAICE_deltaTtherm*            growthNeg       = -SEAICE_deltaTtherm*
539       &                       FICE(I,J)*recip_QI*AREANm1(I,J,bi,bj)       &                      FICE(I,J)*recip_QI*AREANm1(I,J,bi,bj)
540  C     negative growth must not exceed effective ice thickness (=volume)  C     negative growth must not exceed effective ice thickness (=volume)
541  C     (that is, cannot melt more than all the ice)  C     (that is, cannot melt more than all the ice)
542            growthHEFF       = -ONE*MIN(HEFF(I,J,1,bi,bj),growthNeg)            growthHEFF      = -ONE*MIN(HEFF(I,J,bi,bj),growthNeg)
543  C     growthHEFF < 0 means melting  C     growthHEFF < 0 means melting
544            HEFF(I,J,1,bi,bj)= HEFF(I,J,1,bi,bj) + growthHEFF            HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + growthHEFF
545  C     add effective growth to fresh water of ice  C     add effective growth to fresh water of ice
546            saltWtrIce(I,J,bi,bj)  = saltWtrIce(I,J,bi,bj) + growthHEFF            saltWtrIce(I,J,bi,bj) = saltWtrIce(I,J,bi,bj) + growthHEFF
547    
548  C     now calculate QNETI under ice (if any) as the difference between  C     now calculate QNETI under ice (if any) as the difference between
549  C     the available "heat flux" growthNeg and the actual growthHEFF;  C     the available "heat flux" growthNeg and the actual growthHEFF;
# Line 597  C     Note that units of YNEG and frWtrI Line 599  C     Note that units of YNEG and frWtrI
599  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
600  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,
601  CADJ &                        key = iicekey, byte = isbyte  CADJ &                        key = iicekey, byte = isbyte
602  CADJ STORE area(:,:,:,bi,bj) = comlev1_bibj,  CADJ STORE area(:,:,bi,bj)  = comlev1_bibj,
603  CADJ &                        key = iicekey, byte = isbyte  CADJ &                        key = iicekey, byte = isbyte
604  #endif  #endif
605          DO J=1,sNy          DO J=1,sNy
# Line 613  CADJ &                        key = iice Line 615  CADJ &                        key = iice
615  cph)  cph)
616    
617  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
618  CADJ STORE area(:,:,:,bi,bj) = comlev1_bibj,  CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
619  CADJ &                        key = iicekey, byte = isbyte  CADJ &                       key = iicekey, byte = isbyte
620  # ifdef SEAICE_SALINITY  # ifdef SEAICE_SALINITY
621  CADJ STORE hsalt(:,:,bi,bj) = comlev1_bibj,  CADJ STORE hsalt(:,:,bi,bj) = comlev1_bibj,
622  CADJ &                        key = iicekey, byte = isbyte  CADJ &                        key = iicekey, byte = isbyte
# Line 695  C saltWtrIce < 0 : m of sea ice that is Line 697  C saltWtrIce < 0 : m of sea ice that is
697               saltFlux(I,J,bi,bj) =               saltFlux(I,J,bi,bj) =
698       &            HEFFM(I,J,bi,bj)*saltWtrIce(I,J,bi,bj)*       &            HEFFM(I,J,bi,bj)*saltWtrIce(I,J,bi,bj)*
699       &            HSALT(I,J,bi,bj)/       &            HSALT(I,J,bi,bj)/
700       &            (HEFF(I,J,1,bi,bj)-saltWtrIce(I,J,bi,bj))/       &            (HEFF(I,J,bi,bj)-saltWtrIce(I,J,bi,bj))/
701       &            SEAICE_deltaTtherm       &            SEAICE_deltaTtherm
702  #ifdef ALLOW_SALT_PLUME  #ifdef ALLOW_SALT_PLUME
703               saltPlumeFlux(i,j,bi,bj) = 0.0 _d 0               saltPlumeFlux(i,j,bi,bj) = 0.0 _d 0
# Line 707  C update HSALT based on surface saltFlux Line 709  C update HSALT based on surface saltFlux
709            saltFlux(I,J,bi,bj) =            saltFlux(I,J,bi,bj) =
710       &         saltFlux(I,J,bi,bj) + saltFluxAdjust(I,J)       &         saltFlux(I,J,bi,bj) + saltFluxAdjust(I,J)
711  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
712            IF ( HEFF(I,J,1,bi,bj) .EQ. 0.0 ) THEN            IF ( HEFF(I,J,bi,bj) .EQ. 0.0 ) THEN
713               saltFlux(I,J,bi,bj) = saltFlux(I,J,bi,bj) -               saltFlux(I,J,bi,bj) = saltFlux(I,J,bi,bj) -
714       &            HEFFM(I,J,bi,bj) * HSALT(I,J,bi,bj) /       &            HEFFM(I,J,bi,bj) * HSALT(I,J,bi,bj) /
715       &            SEAICE_deltaTtherm       &            SEAICE_deltaTtherm
# Line 736  C NOW GET TOTAL QNET AND QSW Line 738  C NOW GET TOTAL QNET AND QSW
738           IF ( DIAGNOSTICS_IS_ON('SIqneto ',myThid) ) THEN           IF ( DIAGNOSTICS_IS_ON('SIqneto ',myThid) ) THEN
739            DO J=1,sNy            DO J=1,sNy
740             DO I=1,sNx             DO I=1,sNx
741              DIAGarray(I,J) = QNETO(I,J) * (ONE-AREA(I,J,2,bi,bj))              DIAGarray(I,J) = QNETO(I,J) * (ONE-areaNm1(I,J,bi,bj))
742             ENDDO             ENDDO
743            ENDDO            ENDDO
744            CALL DIAGNOSTICS_FILL(DIAGarray,'SIqneto ',0,1,3,bi,bj,myThid)            CALL DIAGNOSTICS_FILL(DIAGarray,'SIqneto ',0,1,3,bi,bj,myThid)
# Line 744  C NOW GET TOTAL QNET AND QSW Line 746  C NOW GET TOTAL QNET AND QSW
746           IF ( DIAGNOSTICS_IS_ON('SIqneti ',myThid) ) THEN           IF ( DIAGNOSTICS_IS_ON('SIqneti ',myThid) ) THEN
747            DO J=1,sNy            DO J=1,sNy
748             DO I=1,sNx             DO I=1,sNx
749              DIAGarray(I,J) = QNETI(I,J) * AREA(I,J,2,bi,bj)              DIAGarray(I,J) = QNETI(I,J) * areaNm1(I,J,bi,bj)
750             ENDDO             ENDDO
751            ENDDO            ENDDO
752            CALL DIAGNOSTICS_FILL(DIAGarray,'SIqneti ',0,1,3,bi,bj,myThid)            CALL DIAGNOSTICS_FILL(DIAGarray,'SIqneti ',0,1,3,bi,bj,myThid)
# Line 789  crg Added by Ralf Giering: do we need DO Line 791  crg Added by Ralf Giering: do we need DO
791  #define DO_WE_NEED_THIS  #define DO_WE_NEED_THIS
792  C NOW ZERO OUTSIDE POINTS  C NOW ZERO OUTSIDE POINTS
793  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
794  CADJ STORE area(:,:,:,bi,bj) = comlev1_bibj,  CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
795  CADJ &                         key = iicekey, byte = isbyte  CADJ &                       key = iicekey, byte = isbyte
796  CADJ STORE heff(:,:,:,bi,bj) = comlev1_bibj,  CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,
797  CADJ &                         key = iicekey, byte = isbyte  CADJ &                       key = iicekey, byte = isbyte
798  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
799          DO J=1,sNy          DO J=1,sNy
800           DO I=1,sNx           DO I=1,sNx
801  C NOW SET AREA(I,J,1,bi,bj)=0 WHERE NO ICE IS  C NOW SET AREA(I,J,bi,bj)=0 WHERE NO ICE IS
802            AREA(I,J,1,bi,bj)=MIN(AREA(I,J,1,bi,bj)            AREA(I,J,bi,bj)=MIN(AREA(I,J,bi,bj)
803       &                         ,HEFF(I,J,1,bi,bj)/.0001 _d 0)       &                         ,HEFF(I,J,bi,bj)/.0001 _d 0)
804           ENDDO           ENDDO
805          ENDDO          ENDDO
806  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
807  CADJ STORE area(:,:,:,bi,bj) = comlev1_bibj,  CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
808  CADJ &                         key = iicekey, byte = isbyte  CADJ &                       key = iicekey, byte = isbyte
809  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
810          DO J=1,sNy          DO J=1,sNy
811           DO I=1,sNx           DO I=1,sNx
812  C NOW TRUNCATE AREA  C NOW TRUNCATE AREA
813  #ifdef DO_WE_NEED_THIS  #ifdef DO_WE_NEED_THIS
814            AREA(I,J,1,bi,bj)=MIN(ONE,AREA(I,J,1,bi,bj))            AREA(I,J,bi,bj)=MIN(ONE,AREA(I,J,bi,bj))
815           ENDDO           ENDDO
816          ENDDO          ENDDO
817  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
818  CADJ STORE area(:,:,:,bi,bj) = comlev1_bibj,  CADJ STORE area(:,:,bi,bj)  = comlev1_bibj,
819  CADJ &                         key = iicekey, byte = isbyte  CADJ &                        key = iicekey, byte = isbyte
820  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,
821  CADJ &                         key = iicekey, byte = isbyte  CADJ &                         key = iicekey, byte = isbyte
822  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
823          DO J=1,sNy          DO J=1,sNy
824           DO I=1,sNx           DO I=1,sNx
825            AREA(I,J,1,bi,bj) = MAX(ZERO,AREA(I,J,1,bi,bj))            AREA(I,J,bi,bj) = MAX(ZERO,AREA(I,J,bi,bj))
826            HSNOW(I,J,bi,bj)  = MAX(ZERO,HSNOW(I,J,bi,bj))            HSNOW(I,J,bi,bj)  = MAX(ZERO,HSNOW(I,J,bi,bj))
827  #endif /* DO_WE_NEED_THIS */  #endif /* DO_WE_NEED_THIS */
828            AREA(I,J,1,bi,bj) = AREA(I,J,1,bi,bj)*HEFFM(I,J,bi,bj)            AREA(I,J,bi,bj) = AREA(I,J,bi,bj)*HEFFM(I,J,bi,bj)
829            HEFF(I,J,1,bi,bj) = HEFF(I,J,1,bi,bj)*HEFFM(I,J,bi,bj)            HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj)*HEFFM(I,J,bi,bj)
830  #ifdef SEAICE_CAP_HEFF  #ifdef SEAICE_CAP_HEFF
831            HEFF(I,J,1,bi,bj)=MIN(MAX_HEFF,HEFF(I,J,1,bi,bj))            HEFF(I,J,bi,bj)=MIN(MAX_HEFF,HEFF(I,J,bi,bj))
832  #endif /* SEAICE_CAP_HEFF */  #endif /* SEAICE_CAP_HEFF */
833            HSNOW(I,J,bi,bj)  = HSNOW(I,J,bi,bj)*HEFFM(I,J,bi,bj)            HSNOW(I,J,bi,bj)  = HSNOW(I,J,bi,bj)*HEFFM(I,J,bi,bj)
834           ENDDO           ENDDO
# Line 838  CADJ &                         key = iic Line 840  CADJ &                         key = iic
840  C     use (abuse) GHEFF to diagnose the total thermodynamic growth rate  C     use (abuse) GHEFF to diagnose the total thermodynamic growth rate
841            DO J=1,sNy            DO J=1,sNy
842             DO I=1,sNx             DO I=1,sNx
843              GHEFF(I,J) = (HEFF(I,J,1,bi,bj)-HEFFNm1(I,J,bi,bj))              GHEFF(I,J) = (HEFF(I,J,bi,bj)-HEFFNm1(I,J,bi,bj))
844       &           /SEAICE_deltaTtherm       &           /SEAICE_deltaTtherm
845             ENDDO             ENDDO
846            ENDDO            ENDDO
# Line 853  C     convert snow to ice if submerged Line 855  C     convert snow to ice if submerged
855           DO J=1,sNy           DO J=1,sNy
856            DO I=1,sNx            DO I=1,sNx
857             hDraft     = (HSNOW(I,J,bi,bj)*SEAICE_rhoSnow             hDraft     = (HSNOW(I,J,bi,bj)*SEAICE_rhoSnow
858       &              +HEFF(I,J,1,bi,bj)*SEAICE_rhoIce)/1000. _d 0       &              +HEFF(I,J,bi,bj)*SEAICE_rhoIce)/1000. _d 0
859  C     here GEFF is the gain of ice due to flooding  C     here GEFF is the gain of ice due to flooding
860             GHEFF(I,J) = hDraft - MIN(hDraft,HEFF(I,J,1,bi,bj))             GHEFF(I,J) = hDraft - MIN(hDraft,HEFF(I,J,bi,bj))
861             HEFF(I,J,1,bi,bj) = HEFF(I,J,1,bi,bj) + GHEFF(I,J)             HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + GHEFF(I,J)
862             HSNOW(I,J,bi,bj)  = MAX(0. _d 0,             HSNOW(I,J,bi,bj)  = MAX(0. _d 0,
863       &          HSNOW(I,J,bi,bj)-GHEFF(I,J)*ICE2SNOW)       &          HSNOW(I,J,bi,bj)-GHEFF(I,J)*ICE2SNOW)
864            ENDDO            ENDDO
# Line 880  C     turn GHEFF into a rate Line 882  C     turn GHEFF into a rate
882          IF ( useRealFreshWaterFlux ) THEN          IF ( useRealFreshWaterFlux ) THEN
883           DO J=1,sNy           DO J=1,sNy
884            DO I=1,sNx            DO I=1,sNx
885             sIceLoad(i,j,bi,bj) = HEFF(I,J,1,bi,bj)*SEAICE_rhoIce             sIceLoad(i,j,bi,bj) = HEFF(I,J,bi,bj)*SEAICE_rhoIce
886       &                         + HSNOW(I,J,bi,bj)*SEAICE_rhoSnow       &                         + HSNOW(I,J,bi,bj)*SEAICE_rhoSnow
887            ENDDO            ENDDO
888           ENDDO           ENDDO
# Line 890  C     turn GHEFF into a rate Line 892  C     turn GHEFF into a rate
892  C     Sources and sinks for sea ice age  C     Sources and sinks for sea ice age
893          DO J=1,sNy          DO J=1,sNy
894           DO I=1,sNx           DO I=1,sNx
895            IF ( AREA(I,J,1,bi,bj) .GT. 0.15 ) THEN            IF ( AREA(I,J,bi,bj) .GT. 0.15 ) THEN
896             IceAge(i,j,bi,bj) = IceAge(i,j,bi,bj) +             IceAge(i,j,bi,bj) = IceAge(i,j,bi,bj) +
897       &          AREA(I,J,1,bi,bj) * SEAICE_deltaTtherm       &          AREA(I,J,bi,bj) * SEAICE_deltaTtherm
898            ELSE            ELSE
899             IceAge(i,j,bi,bj) = ZERO             IceAge(i,j,bi,bj) = ZERO
900            ENDIF            ENDIF
Legend:
Removed from v.1.60  
changed lines
  Added in v.1.61
  ViewVC Help
Powered by ViewVC 1.1.22