New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
limthd_sal.F90 in branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90 @ 3963

Last change on this file since 3963 was 3963, checked in by clem, 11 years ago

bugs correction + creation of glob_max and glob_min in lib_fortran.F90, see ticket:#1116
  • Property svn:keywords set to Id
File size: 12.4 KB
Line 
1MODULE limthd_sal
2   !!======================================================================
3   !!                       ***  MODULE limthd_sal ***
4   !! LIM-3 sea-ice :  computation of salinity variations in the ice
5   !!======================================================================
6   !! History :   -   ! 2003-05 (M. Vancoppenolle) UCL-ASTR first coding for LIM3-1D
7   !!            3.0  ! 2005-12 (M. Vancoppenolle) adapted to the 3-D version
8   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation
9   !!---------------------------------------------------------------------
10#if defined key_lim3
11   !!----------------------------------------------------------------------
12   !!   'key_lim3'                                      LIM-3 sea-ice model
13   !!----------------------------------------------------------------------
14   !!   lim_thd_sal : salinity variations in the ice
15   !!----------------------------------------------------------------------
16   USE par_oce          ! ocean parameters
17   USE phycst           ! physical constants (ocean directory)
18   USE sbc_oce          ! Surface boundary condition: ocean fields
19   USE ice              ! LIM variables
20   USE par_ice          ! LIM parameters
21   USE thd_ice          ! LIM thermodynamics
22   USE limvar           ! LIM variables
23   USE in_out_manager   ! I/O manager
24   USE lib_mpp          ! MPP library
25   USE wrk_nemo         ! work arrays
26   USE lib_fortran      ! to use key_nosignedzero
27
28   IMPLICIT NONE
29   PRIVATE
30
31   PUBLIC   lim_thd_sal        ! called by limthd module
32   PUBLIC   lim_thd_sal_init   ! called by iceini module
33
34   REAL(wp) ::   epsi20 = 1e-20_wp   ! constant values
35   !!----------------------------------------------------------------------
36   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
37   !! $Id$
38   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
39   !!----------------------------------------------------------------------
40CONTAINS
41
42   SUBROUTINE lim_thd_sal( kideb, kiut )
43      !!-------------------------------------------------------------------
44      !!                ***  ROUTINE lim_thd_sal  ***   
45      !!   
46      !! ** Purpose :   computes new salinities in the ice
47      !!
48      !! ** Method  :  4 possibilities
49      !!               -> num_sal = 1 -> constant salinity for z,t
50      !!               -> num_sal = 2 -> S = S(z,t) [simple Vancoppenolle et al 2005]
51      !!               -> num_sal = 3 -> S = S(z)   [multiyear ice]
52      !!               -> num_sal = 4 -> S = S(h)   [Cox and Weeks 74]
53      !!---------------------------------------------------------------------
54      INTEGER, INTENT(in) ::  kideb, kiut   ! thickness category index
55      !
56      INTEGER  ::   ji, jk     ! dummy loop indices
57      INTEGER  ::   zji, zjj   ! local integers
58      REAL(wp) ::   zsold, iflush, iaccrbo, igravdr, isnowic, i_ice_switch,  ztmelts   ! local scalars
59      REAL(wp) ::   zaaa, zbbb, zccc, zdiscrim   ! local scalars
60      REAL(wp), POINTER, DIMENSION(:) ::   ze_init, zhiold, zsiold
61      !!---------------------------------------------------------------------
62
63      CALL wrk_alloc( jpij, ze_init, zhiold, zsiold )
64
65      !------------------------------------------------------------------------------|
66      ! 1) Constant salinity, constant in time                                       |
67      !------------------------------------------------------------------------------|
68!!gm comment: if num_sal = 1 s_i_b and sm_i_b can be set to bulk_sal one for all in the initialisation phase !!
69      IF( num_sal == 1 ) THEN
70         !
71         s_i_b  (kideb:kiut,1:nlay_i) =  bulk_sal 
72         sm_i_b (kideb:kiut)          =  bulk_sal
73         s_i_new(kideb:kiut)          =  bulk_sal
74        !
75      ENDIF
76
77      !------------------------------------------------------------------------------|
78      !  Module 2 : Constant salinity varying in time                                |
79      !------------------------------------------------------------------------------|
80
81      IF(  num_sal == 2  .OR.  num_sal == 4  ) THEN
82
83         !---------------------------------
84         ! Thickness at previous time step
85         !---------------------------------
86         DO ji = kideb, kiut
87            zhiold(ji) = ht_i_b(ji) - dh_i_bott(ji) - dh_snowice(ji) - dh_i_surf(ji)
88            zsiold(ji) = sm_i_b(ji)
89         END DO
90
91         !---------------------
92         ! Global heat content
93         !---------------------
94         ze_init(:)  =  0._wp
95         DO jk = 1, nlay_i
96            DO ji = kideb, kiut
97               ze_init(ji) = ze_init(ji) + q_i_b(ji,jk) * ht_i_b(ji) / REAL (nlay_i )
98            END DO
99         END DO
100
101         DO ji = kideb, kiut
102            !
103            ! Switches
104            !----------
105            iflush       =         MAX( 0._wp , SIGN( 1.0 , t_su_b(ji) - rtt )        )    ! =1 if summer
106            igravdr      =         MAX( 0._wp , SIGN( 1.0 , t_bo_b(ji) - t_su_b(ji) ) )    ! =1 if t_su < t_bo
107            iaccrbo      =         MAX( 0._wp , SIGN( 1.0 , dh_i_bott(ji) )           )    ! =1 if bottom accretion
108            i_ice_switch = 1._wp - MAX ( 0._wp , SIGN( 1._wp , - ht_i_b(ji) + 1.e-2 ) )
109            isnowic      = 1._wp - MAX ( 0._wp , SIGN( 1._wp , - dh_snowice(ji) ) ) * i_ice_switch   ! =1 if snow ice formation
110
111            !---------------------
112            ! Salinity tendencies
113            !---------------------
114            !                                   ! drainage by gravity drainage
115            dsm_i_gd_1d(ji) = - igravdr * MAX( sm_i_b(ji) - sal_G , 0._wp ) / time_G * rdt_ice 
116            !                                   ! drainage by flushing 
117            dsm_i_fl_1d(ji) = - iflush * MAX( sm_i_b(ji) - sal_F , 0._wp ) / time_F * rdt_ice
118
119            !-----------------
120            ! Update salinity   
121            !-----------------
122            ! only drainage terms ( gravity drainage and flushing )
123            ! snow ice / bottom sources are added in lim_thd_ent to conserve energy
124            sm_i_b(ji) = sm_i_b(ji) + dsm_i_fl_1d(ji) + dsm_i_gd_1d(ji)
125
126            ! if no ice, salinity = 0.1
127            i_ice_switch = 1._wp - MAX ( 0._wp, SIGN( 1._wp , - ht_i_b(ji) ) )
128            sm_i_b(ji)   = i_ice_switch * sm_i_b(ji) + s_i_min * ( 1._wp - i_ice_switch )
129
130            !----------------------------
131            ! Heat flux - brine drainage
132            !----------------------------
133            fhbri_1d(ji) = 0._wp
134
135            !----------------------------
136            ! Salt flux - brine drainage
137            !----------------------------
138            fsbri_1d(ji) = fsbri_1d(ji) - i_ice_switch * rhoic * a_i_b(ji) * ht_i_b(ji) * ( sm_i_b(ji) - zsiold(ji) ) / rdt_ice
139            IF( num_sal == 4 ) fsbri_1d(ji) = 0._wp
140
141        END DO ! ji
142
143         ! Salinity profile
144         CALL lim_var_salprof1d( kideb, kiut )
145
146         !----------------------------
147         ! Heat flux - brine drainage
148         !----------------------------
149
150!clem:move         DO ji = kideb, kiut
151!!gm useless
152            ! iflush  : 1 if summer
153!clem            iflush  =  MAX( 0._wp , SIGN ( 1._wp , t_su_b(ji) - rtt ) )
154            ! igravdr : 1 if t_su lt t_bo
155!clem            igravdr =  MAX( 0._wp , SIGN ( 1._wp , t_bo_b(ji) - t_su_b(ji) ) )
156            ! iaccrbo : 1 if bottom accretion
157!clem            iaccrbo =  MAX( 0._wp , SIGN ( 1._wp , dh_i_bott(ji) ) )
158!!gm end useless
159            !
160!clem:move            fhbri_1d(ji) = 0._wp
161!clem:move         END DO ! ji
162
163         !----------------------------
164         ! Salt flux - brine drainage
165         !----------------------------
166!clem:move          DO ji = kideb, kiut
167!clem:move            i_ice_switch = 1._wp - MAX ( 0._wp, SIGN( 1._wp , - ht_i_b(ji) ) )
168!clem:move            fsbri_1d(ji) = fsbri_1d(ji) - i_ice_switch * rhoic * a_i_b(ji) * ht_i_b(ji) * ( sm_i_b(ji) - zsiold(ji) ) / rdt_ice
169            !i_ice_switch = 1._wp - MAX ( 0._wp, SIGN( 1._wp , - zhiold(ji) ) )
170            !fsbri_1d(ji) = fsbri_1d(ji) - i_ice_switch * rhoic * a_i_b(ji) * zhiold(ji) * ( sm_i_b(ji) - zsiold(ji) ) / rdt_ice
171            !clem fsbri_1d(ji) = fsbri_1d(ji) - i_ice_switch * rhoic * a_i_b(ji) * ht_i_b(ji)         &
172            !clem     &         * ( MAX(dsm_i_gd_1d(ji) + dsm_i_fl_1d(ji), sm_i_b(ji) - zsiold(ji) ) ) / rdt_ice
173!clem:move            IF( num_sal == 4 ) fsbri_1d(ji) = 0._wp
174!clem:move          END DO ! ji
175
176         ! Only necessary for conservation check since salinity is modified
177         !--------------------
178         ! Temperature update
179         !--------------------
180         DO jk = 1, nlay_i
181            DO ji = kideb, kiut
182               ztmelts    =  -tmut*s_i_b(ji,jk) + rtt
183               !Conversion q(S,T) -> T (second order equation)
184               zaaa         =  cpic
185               zbbb         =  ( rcp - cpic ) * ( ztmelts - rtt ) + q_i_b(ji,jk) / rhoic - lfus
186               zccc         =  lfus * ( ztmelts - rtt )
187               zdiscrim     =  SQRT(  MAX( zbbb*zbbb - 4.0*zaaa*zccc, 0._wp )  )
188               t_i_b(ji,jk) =  rtt - ( zbbb + zdiscrim ) / ( 2.0 *zaaa )
189            END DO
190         END DO
191         !
192      ENDIF ! num_sal .EQ. 2
193
194      !------------------------------------------------------------------------------|
195      !  Module 3 : Profile of salinity, constant in time                            |
196      !------------------------------------------------------------------------------|
197
198      IF( num_sal == 3 )   CALL lim_var_salprof1d( kideb, kiut )
199
200      !------------------------------------------------------------------------------|
201      !  Module 4 : Constant salinity varying in time                                |
202      !------------------------------------------------------------------------------|
203
204      IF( num_sal == 5 ) THEN      ! Cox and Weeks, 1974
205         !
206         DO ji = kideb, kiut
207            zsold = sm_i_b(ji)
208            IF( ht_i_b(ji) < 0.4 ) THEN
209               sm_i_b(ji) = 14.24 - 19.39 * ht_i_b(ji) 
210            ELSE
211               sm_i_b(ji) =  7.88 - 1.59 * ht_i_b(ji)
212               sm_i_b(ji) = MIN( sm_i_b(ji) , zsold ) 
213            ENDIF
214            IF( ht_i_b(ji) > 3.06918239 ) THEN
215               sm_i_b(ji) = 3._wp
216            ENDIF
217            DO jk = 1, nlay_i
218               s_i_b(ji,jk)   = sm_i_b(ji)
219            END DO
220         END DO
221         !
222      ENDIF ! num_sal
223
224      !
225      CALL wrk_dealloc( jpij, ze_init, zhiold, zsiold )
226      !
227   END SUBROUTINE lim_thd_sal
228
229
230   SUBROUTINE lim_thd_sal_init
231      !!-------------------------------------------------------------------
232      !!                  ***  ROUTINE lim_thd_sal_init  ***
233      !!
234      !! ** Purpose :   initialization of ice salinity parameters
235      !!
236      !! ** Method  :   Read the namicesal namelist and check the parameter
237      !!              values called at the first timestep (nit000)
238      !!
239      !! ** input   :   Namelist namicesal
240      !!-------------------------------------------------------------------
241      NAMELIST/namicesal/ num_sal, bulk_sal, sal_G, time_G, sal_F, time_F,   &
242         &                s_i_max, s_i_min, s_i_0, s_i_1
243      !!-------------------------------------------------------------------
244      !
245      REWIND( numnam_ice )                   ! Read Namelist namicesal
246      READ  ( numnam_ice  , namicesal )
247      !
248      IF(lwp) THEN                           ! control print
249         WRITE(numout,*)
250         WRITE(numout,*) 'lim_thd_sal_init : Ice parameters for salinity '
251         WRITE(numout,*) '~~~~~~~~~~~~~~~~'
252         WRITE(numout,*) ' switch for salinity num_sal        : ', num_sal
253         WRITE(numout,*) ' bulk salinity value if num_sal = 1 : ', bulk_sal
254         WRITE(numout,*) ' restoring salinity for GD          : ', sal_G
255         WRITE(numout,*) ' restoring time for GD              : ', time_G
256         WRITE(numout,*) ' restoring salinity for flushing    : ', sal_F
257         WRITE(numout,*) ' restoring time for flushing        : ', time_F
258         WRITE(numout,*) ' Maximum tolerated ice salinity     : ', s_i_max
259         WRITE(numout,*) ' Minimum tolerated ice salinity     : ', s_i_min
260         WRITE(numout,*) ' 1st salinity for salinity profile  : ', s_i_0
261         WRITE(numout,*) ' 2nd salinity for salinity profile  : ', s_i_1
262      ENDIF
263      !
264   END SUBROUTINE lim_thd_sal_init
265
266#else
267   !!----------------------------------------------------------------------
268   !!   Default option         Dummy Module          No LIM-3 sea-ice model
269   !!----------------------------------------------------------------------
270#endif
271   !!======================================================================
272END MODULE limthd_sal
Note: See TracBrowser for help on using the repository browser.