Changeset 4269 for branches/2013/dev_LOCEAN_CMCC_INGV_MERC_UKMO_2013/NEMOGCM/NEMO/LIM_SRC_3

Timestamp:

2013-11-20T09:48:15+01:00 (10 years ago)

Author:

clem

Message:

update lim3 for bdy purpose

Location:

branches/2013/dev_LOCEAN_CMCC_INGV_MERC_UKMO_2013/NEMOGCM/NEMO/LIM_SRC_3

Files:

• limcat_1D.F90 (modified) (5 diffs)
• limrhg.F90 (modified) (2 diffs)

branches/2013/dev_LOCEAN_CMCC_INGV_MERC_UKMO_2013/NEMOGCM/NEMO/LIM_SRC_3/limcat_1D.F90

@@ -45 +45 @@
       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zht_i, zht_s, za_i ! output ice/snow variables
       REAL(wp) :: epsi06 = 1.0e-6
-  
+      REAL(wp) ::   zc1, zc2, zc3, zx1   ! local scalars
+      REAL(wp), DIMENSION(0:jpl)   ::   zhi_max         !:Boundary of ice thickness categories in thickness space
+ 
       IF( nn_timing == 1 )  CALL timing_start('limcat_1D')
       !--------------------------------------------------------------------
@@ -66 +68 @@
       !         fulfills ice volume concervation between input and output (ztests=4) 
       !--------------------------------------------------------------------------------------
+
+      !- Thickness categories boundaries
+      !    hi_max is calculated in iceini.F90 but since limcat_1D.F90 routine 
+      !    is called before (in bdydta.F90), one must recalculate it.   
+      !    Note clem: there may be a way of doing things cleaner
+      !----------------------------------
+      zhi_max(:) = 0._wp
+      zc1 =  3._wp / REAL( jpl , wp ) ; zc2 = 10._wp * zc1 ; zc3 =  3._wp
+      DO jl = 1, jpl
+         zx1 = REAL( jl-1 , wp ) / REAL( jpl , wp )
+         zhi_max(jl) = zhi_max(jl-1) + zc1 + zc2 * ( 1._wp + TANH( zc3 * ( zx1 - 1._wp ) ) )
+      END DO
  
       ! ----------------------------------------
@@ -100 +114 @@
                ! Fill ice thicknesses except the last one (i_fill) by (hmax-hmin)/2 
                DO jl = 1, i_fill - 1
-                  zht_i(ji,jl) = ( hi_max(jl) + hi_max(jl-1) ) / 2.
+                  zht_i(ji,jl) = ( zhi_max(jl) + zhi_max(jl-1) ) / 2.
                END DO
                
@@ -106 +120 @@
                jl0 = i_fill
                DO jl = 1, i_fill
-                  IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN
+                  IF ( ( zhti(ji) >= zhi_max(jl-1) ) .AND. ( zhti(ji) < zhi_max(jl) ) ) THEN
                      jl0 = jl
            CYCLE
@@ -141 +155 @@
             
             ! Test 3: thickness of the last category is in-bounds ?
-            IF ( zht_i(ji,i_fill) >= hi_max(i_fill-1) ) ztest_3 = 1
+            IF ( zht_i(ji,i_fill) >= zhi_max(i_fill-1) ) ztest_3 = 1
             
             ! Test 4: positivity of ice concentrations

branches/2013/dev_LOCEAN_CMCC_INGV_MERC_UKMO_2013/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90

@@ -426 +426 @@
                !-Calculate stress tensor components zs1 and zs2 
                !-at centre of grid cells (see section 3.5 of CICE user's guide).
-               zs1(ji,jj) = ( zs1(ji,jj) - dtotel*( ( 1._wp - alphaevp) * zs1(ji,jj) +   &
-                  &          ( delta / deltat(ji,jj) - zdd(ji,jj) / deltat(ji,jj) ) * zpresh(ji,jj) ) )  &       
-                  &          / ( 1._wp + alphaevp * dtotel )
-
-               zs2(ji,jj) = ( zs2(ji,jj) - dtotel * ( ( 1._wp - alphaevp ) * ecc2 * zs2(ji,jj) -   &
-                             zdt(ji,jj) / deltat(ji,jj) * zpresh(ji,jj) ) )   &
-                  &          / ( 1._wp + alphaevp * ecc2 * dtotel )
+               !zs1(ji,jj) = ( zs1(ji,jj) - dtotel*( ( 1._wp - alphaevp) * zs1(ji,jj) +   &
+               !   &          ( delta / deltat(ji,jj) - zdd(ji,jj) / deltat(ji,jj) ) * zpresh(ji,jj) ) )  &       
+               !   &          / ( 1._wp + alphaevp * dtotel )
+
+               !zs2(ji,jj) = ( zs2(ji,jj) - dtotel * ( ( 1._wp - alphaevp ) * ecc2 * zs2(ji,jj) -   &
+               !              zdt(ji,jj) / deltat(ji,jj) * zpresh(ji,jj) ) )   &
+               !   &          / ( 1._wp + alphaevp * ecc2 * dtotel )
 
                ! new formulation from S. Bouillon to help stabilizing the code (no need of alphaevp)
-               !zs1(ji,jj) = ( zs1(ji,jj) + dtotel * ( ( zdd(ji,jj) / deltat(ji,jj) - delta / deltat(ji,jj) )  &
-               !   &         * zpresh(ji,jj) ) ) / ( 1._wp + dtotel )
-               !zs2(ji,jj) = ( zs2(ji,jj) + dtotel * ( ecci * zdt(ji,jj) / deltat(ji,jj) * zpresh(ji,jj) ) )  &
-               !   &         / ( 1._wp + dtotel )
+               zs1(ji,jj) = ( zs1(ji,jj) + dtotel * ( ( zdd(ji,jj) / deltat(ji,jj) - delta / deltat(ji,jj) )  &
+                  &         * zpresh(ji,jj) ) ) / ( 1._wp + dtotel )
+               zs2(ji,jj) = ( zs2(ji,jj) + dtotel * ( ecci * zdt(ji,jj) / deltat(ji,jj) * zpresh(ji,jj) ) )  &
+                  &         / ( 1._wp + dtotel )
 
             END DO
@@ -472 +472 @@
 
                !-Calculate stress tensor component zs12 at corners (see section 3.5 of CICE user's guide).
-               zs12(ji,jj) = ( zs12(ji,jj) - dtotel * ( (1.0-alphaevp) * ecc2 * zs12(ji,jj) - zds(ji,jj) /  &
-                  &          ( 2._wp * deltac(ji,jj) ) * zpreshc(ji,jj) ) )  &
-                  &          / ( 1._wp + alphaevp * ecc2 * dtotel ) 
+               !zs12(ji,jj) = ( zs12(ji,jj) - dtotel * ( (1.0-alphaevp) * ecc2 * zs12(ji,jj) - zds(ji,jj) /  &
+               !   &          ( 2._wp * deltac(ji,jj) ) * zpreshc(ji,jj) ) )  &
+               !   &          / ( 1._wp + alphaevp * ecc2 * dtotel ) 
 
                ! new formulation from S. Bouillon to help stabilizing the code (no need of alphaevp)
-               !zs12(ji,jj) = ( zs12(ji,jj) + dtotel *  &
-               !   &          ( ecci * zds(ji,jj) / ( 2._wp * deltac(ji,jj) ) * zpreshc(ji,jj) ) )  &
-               !   &          / ( 1.0 + dtotel ) 
+               zs12(ji,jj) = ( zs12(ji,jj) + dtotel *  &
+                  &          ( ecci * zds(ji,jj) / ( 2._wp * deltac(ji,jj) ) * zpreshc(ji,jj) ) )  &
+                  &          / ( 1.0 + dtotel ) 
 
             END DO ! ji