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 trunk/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90 @ 2528

Last change on this file since 2528 was 2528, checked in by rblod, 13 years ago

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