Changeset 4833 for trunk/NEMOGCM

Timestamp:

2014-11-06T15:58:37+01:00 (9 years ago)

Author:

flavoni

Message:

fix LIM3 reproducibility problems, see ticket #1398

File:

• trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90 (modified) (12 diffs)

trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90

@@ -29 +29 @@
    USE lib_mpp        ! MPP library
    USE wrk_nemo       ! work arrays
+   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
    USE limthd_ent
@@ -133 +134 @@
                   !Energy of melting q(S,T) [J.m-3]
                   zindb = 1._wp - MAX(  0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) + epsi10 )  )   !0 if no ice and 1 if yes
-                  e_i(ji,jj,jk,jl) = zindb * e_i(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_i(ji,jj,jl) ,  epsi10 ) ) * REAL( nlay_i )
+                  e_i(ji,jj,jk,jl) = zindb * e_i(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_i(ji,jj,jl) ,  epsi10 ) ) * REAL( nlay_i, wp )
                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac
                END DO
@@ -171 +172 @@
          zgamafr = 0.03
 
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-
+         DO jj = 2, jpj
+            DO ji = 2, jpi
                IF ( qlead(ji,jj) < 0._wp ) THEN
                   !-------------
@@ -243 +243 @@
             END DO ! loop on ji ends
          END DO ! loop on jj ends
+      ! 
+      CALL lbc_lnk( zvrel(:,:), 'T', 1. )
+      CALL lbc_lnk( hicol(:,:), 'T', 1. )
 
       ENDIF ! End of computation of frazil ice collection thickness
@@ -255 +258 @@
       ! This occurs if open water energy budget is negative
       nbpac = 0
+      npac(:) = 0
+      !
       DO jj = 1, jpj
          DO ji = 1, jpi
@@ -315 +320 @@
          ! Keep old ice areas and volume in memory
          !-----------------------------------------
-         zv_old(:,:) = zv_i_1d(:,:) 
-         za_old(:,:) = za_i_1d(:,:)
-
+         zv_old(1:nbpac,:) = zv_i_1d(1:nbpac,:) 
+         za_old(1:nbpac,:) = za_i_1d(1:nbpac,:)
          !----------------------
          ! Thickness of new ice
@@ -324 +328 @@
             zh_newice(ji) = hiccrit
          END DO
-         IF( fraz_swi == 1 ) zh_newice(:) = hicol_b(:)
+         IF( fraz_swi == 1 ) zh_newice(1:nbpac) = hicol_b(1:nbpac)
 
          !----------------------
@@ -331 +335 @@
          SELECT CASE ( num_sal )
          CASE ( 1 )                    ! Sice = constant 
-            zs_newice(:) = bulk_sal
+            zs_newice(1:nbpac) = bulk_sal
          CASE ( 2 )                    ! Sice = F(z,t) [Vancoppenolle et al (2005)]
             DO ji = 1, nbpac
@@ -339 +343 @@
             END DO
          CASE ( 3 )                    ! Sice = F(z) [multiyear ice]
-            zs_newice(:) =   2.3
+            zs_newice(1:nbpac) =   2.3
          END SELECT
 
@@ -472 +476 @@
                za_i_1d(ji,jl) = zinda * za_i_1d(ji,jl)               
                zv_i_1d(ji,jl) = zv_i_1d(ji,jl) + zv_newfra
-
                ! for remapping
                h_i_old (ji,nlay_i+1) = zv_newfra
@@ -479 +482 @@
 
             ! --- Ice enthalpy remapping --- !
-            IF( zv_newfra > 0._wp ) THEN
-               CALL lim_thd_ent( 1, nbpac, ze_i_1d(1:nbpac,:,jl) ) 
-            ENDIF
-
+            CALL lim_thd_ent( 1, nbpac, ze_i_1d(1:nbpac,:,jl) ) 
          ENDDO
 
@@ -534 +534 @@
                DO ji = 1, jpi
                   ! heat content in Joules
-                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * area(ji,jj) * v_i(ji,jj,jl) / ( REAL( nlay_i ) * unit_fac ) 
+                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * area(ji,jj) * v_i(ji,jj,jl) / ( REAL( nlay_i ,wp ) * unit_fac ) 
                END DO
             END DO