Changeset 4332 for branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limcat_1D.F90

Timestamp:

2013-12-11T15:38:42+01:00 (10 years ago)

Author:

clem

Message:

update LIM3 to fix remaining bugs. Now working in global and regional config.

File:

• branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limcat_1D.F90 (modified) (11 diffs)

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

@@ -44 +44 @@
       REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zai    ! input ice/snow variables
       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zht_i, zht_s, za_i ! output ice/snow variables
-      REAL(wp) :: epsi06 = 1.0e-6
-  
+      REAL(wp)                     ::   epsi06 = 1.0e-6
+      REAL(wp)                     ::   zc1, zc2, zc3, zx1, zdh   ! 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')
       !--------------------------------------------------------------------
@@ -51 +53 @@
       !--------------------------------------------------------------------
       ijpij = SIZE(zhti,1)
-      zht_i(1:ijpij,1:jpl) = 0.d0
-      zht_s(1:ijpij,1:jpl) = 0.d0
-      za_i (1:ijpij,1:jpl) = 0.d0
+      zht_i(1:ijpij,1:jpl) = 0._wp
+      zht_s(1:ijpij,1:jpl) = 0._wp
+      za_i (1:ijpij,1:jpl) = 0._wp
 
       !------------------------------------------------------------------------------------
@@ -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
  
       ! ----------------------------------------
@@ -71 +85 @@
       ! ----------------------------------------
       DO ji = 1, ijpij
-        ! snow thickness in each category 
-         zht_s(ji,1:jpl) = zhts(ji)
          
+         IF( zhti(ji) > 0._wp ) THEN
+
          ! initialisation of tests
          ztest_1 = 0
@@ -87 +101 @@
             
             ! initialisation of ice variables for each try
-            zht_i(ji,1:jpl) = 0.d0
-            za_i (ji,1:jpl) = 0.d0
+            zht_i(ji,1:jpl) = 0._wp
+            za_i (ji,1:jpl) = 0._wp
             
             ! *** case very thin ice: fill only category 1
@@ -94 +108 @@
                zht_i(ji,1) = zhti(ji)
                za_i (ji,1) = zai (ji)
-               
                ! *** case ice is thicker: fill categories >1
             ELSE
@@ -100 +113 @@
                ! 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) ) * 0.5_wp
                END DO
                
@@ -106 +119 @@
                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
@@ -116 +129 @@
                DO jl = 1, i_fill - 1
                   IF ( jl == jl0 ) CYCLE
-                  zarg           = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) / 2. )
+                  zarg           = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )
                   za_i(ji,jl) =   za_i (ji,jl0) * EXP(-zarg**2)
                END DO
@@ -141 +154 @@
             
             ! 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
             ztest_4 = 1
             DO jl = 1, i_fill
-               IF ( za_i(ji,jl) < 0.0d0 ) ztest_4 = 0
+               IF ( za_i(ji,jl) < 0._wp ) ztest_4 = 0
             END DO
             
@@ -164 +177 @@
          !   WRITE(numout,*) 'za_i(ji,jpl)=',za_i(ji,:)
          !ENDIF
-            
+         
+         ENDIF ! if zhti > 0
+
       END DO ! i loop
+
+      ! ------------------------------------------------
+      ! Adding Snow in each category where za_i is not 0
+      ! ------------------------------------------------ 
+      DO jl = 1, jpl
+         DO ji = 1, ijpij
+            IF( za_i(ji,jl) > 0._wp ) THEN
+               zht_s(ji,jl) = zht_i(ji,jl) * ( zhts(ji) / zhti(ji) )
+               ! In case snow load is in excess that would lead to transformation from snow to ice
+               ! Then, transfer the snow excess into the ice (different from limthd_dh)
+               zdh = MAX( 0._wp, ( rhosn * zht_s(ji,jl) + ( rhoic - rau0 ) * zht_i(ji,jl) ) * r1_rau0 ) 
+               ! recompute ht_i, ht_s avoiding out of bounds values
+               zht_i(ji,jl) = MIN( zhi_max(jl), zht_i(ji,jl) + zdh )
+               zht_s(ji,jl) = MAX( 0._wp, zht_s(ji,jl) - zdh * rhoic / rhosn )
+            ENDIF
+         ENDDO
+      ENDDO
 
       IF( nn_timing == 1 )  CALL timing_stop('limcat_1D')