Changeset 8325 for branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90

Timestamp:

2017-07-12T16:51:20+02:00 (7 years ago)

Author:

clem

Message:

STEP4 (1): put all thermodynamics in 1D

File:

• branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90 (modified) (7 diffs)

branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90

@@ -131 +131 @@
       CALL lbc_lnk( zfric, 'T',  1. )
       !
-      !----------------------------------!
-      ! Initialization and units change
-      !----------------------------------!
-      ftr_ice(:,:,:) = 0._wp  ! part of solar radiation transmitted through the ice
-
-      ! Change the units of heat content; from J/m2 to J/m3
-      DO jl = 1, jpl
-         DO jk = 1, nlay_i
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  rswitch = MAX(  0._wp , SIGN( 1._wp , v_i(ji,jj,jl) - epsi20 )  )
-                  !Energy of melting q(S,T) [J.m-3]
-                  e_i(ji,jj,jk,jl) = rswitch * e_i(ji,jj,jk,jl) / MAX( v_i(ji,jj,jl) , epsi20 ) * REAL( nlay_i )
-               END DO
-            END DO
-         END DO
-         DO jk = 1, nlay_s
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  rswitch = MAX(  0._wp , SIGN( 1._wp , v_s(ji,jj,jl) - epsi20 )  )
-                  !Energy of melting q(S,T) [J.m-3]
-                  e_s(ji,jj,jk,jl) = rswitch * e_s(ji,jj,jk,jl) / MAX( v_s(ji,jj,jl) , epsi20 ) * REAL( nlay_s )
-               END DO
-            END DO
-         END DO
-      END DO
+      ftr_ice(:,:,:) = 0._wp  ! initialization (part of solar radiation transmitted through the ice)
 
       !--------------------------------------------------------------------!
@@ -269 +244 @@
                               CALL lim_thd_1d2d( nbpb, jl, 1 )               ! --- Move to 1D arrays --- !
             !
+            DO jk = 1, nlay_i                                                ! --- Change units from J/m2 to J/m3 --- !
+               WHERE( ht_i_1d(:)>0._wp ) e_i_1d(:,jk) = e_i_1d(:,jk) / (ht_i_1d(:) * a_i_1d(:)) * nlay_i
+            ENDDO
+            DO jk = 1, nlay_s
+               WHERE( ht_s_1d(:)>0._wp ) e_s_1d(:,jk) = e_s_1d(:,jk) / (ht_s_1d(:) * a_i_1d(:)) * nlay_s
+            ENDDO
+            !
             IF( ln_limdH )    CALL lim_thd_dif( 1, nbpb )                    ! --- Ice/Snow Temperature profile --- !
             !
             IF( ln_limdH )    CALL lim_thd_dh( 1, nbpb )                     ! --- Ice/Snow thickness --- !    
             !
-            IF( ln_limdH )    CALL lim_thd_ent( 1, nbpb, q_i_1d(1:nbpb,:) )  ! --- Ice enthalpy remapping --- !
+            IF( ln_limdH )    CALL lim_thd_ent( 1, nbpb, e_i_1d(1:nbpb,:) )  ! --- Ice enthalpy remapping --- !
             !
                               CALL lim_thd_sal( 1, nbpb )                    ! --- Ice salinity --- !    
@@ -285 +267 @@
             END IF
             !
+            DO jk = 1, nlay_i                                                ! --- Change units from J/m3 to J/m2 --- !
+               e_i_1d(:,jk) = e_i_1d(:,jk) * ht_i_1d(:) * a_i_1d(:) * r1_nlay_i
+            ENDDO
+            DO jk = 1, nlay_s
+               e_s_1d(:,jk) = e_s_1d(:,jk) * ht_s_1d(:) * a_i_1d(:) * r1_nlay_s
+            ENDDO
+            !
                               CALL lim_thd_1d2d( nbpb, jl, 2 )               ! --- Move to 2D arrays --- !
             !
@@ -294 +283 @@
       IF( ln_limdA)           CALL lim_thd_da                                ! --- lateral melting --- !
 
-      ! Enthalpies are global variables we have to readjust the units (heat content in J/m2)
-      DO jl = 1, jpl
-         DO jk = 1, nlay_i
-            e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * a_i(:,:,jl) * ht_i(:,:,jl) * r1_nlay_i
-         END DO
-         DO jk = 1, nlay_s
-            e_s(:,:,jk,jl) = e_s(:,:,jk,jl) * a_i(:,:,jl) * ht_s(:,:,jl) * r1_nlay_s
-         END DO
-      END DO
- 
       ! Change thickness to volume
       v_i(:,:,:)   = ht_i(:,:,:) * a_i(:,:,:)
@@ -375 +354 @@
             ! Conversion q(S,T) -> T (second order equation)
             zaaa          =  cpic
-            zbbb          =  ( rcp - cpic ) * ( ztmelts - rt0 ) + q_i_1d(ji,jk) * r1_rhoic - lfus
+            zbbb          =  ( rcp - cpic ) * ( ztmelts - rt0 ) + e_i_1d(ji,jk) * r1_rhoic - lfus
             zccc          =  lfus * ( ztmelts - rt0 )
             zdiscrim      =  SQRT( MAX( zbbb * zbbb - 4._wp * zaaa * zccc, 0._wp ) )
@@ -451 +430 @@
          DO jk = 1, nlay_s
             CALL tab_2d_1d( nbpb, t_s_1d(1:nbpb,jk), t_s(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) )
-            CALL tab_2d_1d( nbpb, q_s_1d(1:nbpb,jk), e_s(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) )
+            CALL tab_2d_1d( nbpb, e_s_1d(1:nbpb,jk), e_s(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) )
          END DO
          DO jk = 1, nlay_i
             CALL tab_2d_1d( nbpb, t_i_1d(1:nbpb,jk), t_i(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) )
-            CALL tab_2d_1d( nbpb, q_i_1d(1:nbpb,jk), e_i(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) )
+            CALL tab_2d_1d( nbpb, e_i_1d(1:nbpb,jk), e_i(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) )
             CALL tab_2d_1d( nbpb, s_i_1d(1:nbpb,jk), s_i(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) )
          END DO
@@ -523 +502 @@
          DO jk = 1, nlay_s
             CALL tab_1d_2d( nbpb, t_s(:,:,jk,jl), npb, t_s_1d     (1:nbpb,jk), jpi, jpj)
-            CALL tab_1d_2d( nbpb, e_s(:,:,jk,jl), npb, q_s_1d     (1:nbpb,jk), jpi, jpj)
+            CALL tab_1d_2d( nbpb, e_s(:,:,jk,jl), npb, e_s_1d     (1:nbpb,jk), jpi, jpj)
          END DO
          DO jk = 1, nlay_i
             CALL tab_1d_2d( nbpb, t_i(:,:,jk,jl), npb, t_i_1d     (1:nbpb,jk), jpi, jpj)
-            CALL tab_1d_2d( nbpb, e_i(:,:,jk,jl), npb, q_i_1d     (1:nbpb,jk), jpi, jpj)
+            CALL tab_1d_2d( nbpb, e_i(:,:,jk,jl), npb, e_i_1d     (1:nbpb,jk), jpi, jpj)
             CALL tab_1d_2d( nbpb, s_i(:,:,jk,jl), npb, s_i_1d     (1:nbpb,jk), jpi, jpj)
          END DO