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 tags/nemo_v3_2_2/NEMO/LIM_SRC_3 – NEMO

source: tags/nemo_v3_2_2/NEMO/LIM_SRC_3/limthd_sal.F90 @ 3532

Last change on this file since 3532 was 2477, checked in by cetlod, 13 years ago

v3.2:remove hardcoded value of num_sal in limrst.F90, see ticket #633
  • Property svn:keywords set to Id
File size: 14.6 KB
Line 
1MODULE limthd_sal
2   !!----------------------------------------------------------------------
3   !!   'key_lim3'                                      LIM3 sea-ice model
4   !!----------------------------------------------------------------------
5   !!======================================================================
6   !!                       ***  MODULE limthd_sal ***
7   !!                   computation salinity variations in
8   !!                               the ice
9   !!======================================================================
10#if defined key_lim3
11   !!----------------------------------------------------------------------
12   !!   lim_thd_sal : salinity variations in the ice
13   !! * Modules used
14   USE par_oce          ! ocean parameters
15   USE phycst           ! physical constants (ocean directory)
16   USE sbc_oce          ! Surface boundary condition: ocean fields
17   USE thd_ice
18   USE ice
19   USE in_out_manager
20   USE limvar
21   USE par_ice
22
23   IMPLICIT NONE
24   PRIVATE
25
26   !! * Routine accessibility
27   PUBLIC lim_thd_sal        ! called by lim_thd
28   PUBLIC lim_thd_sal_init   ! called by lim_thd
29
30   !! * Module variables
31   REAL(wp)  ::           &  ! constant values
32      epsi20 = 1e-20   ,  &
33      epsi13 = 1e-13   ,  &
34      zzero  = 0.e0    ,  &
35      zone   = 1.e0
36
37   !!----------------------------------------------------------------------
38   !!   LIM 3.0,  UCL-ASTR-LOCEAN-IPSL (2008)
39   !! $Id$
40   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
41   !!----------------------------------------------------------------------
42
43CONTAINS
44
45   SUBROUTINE lim_thd_sal(kideb,kiut)
46      !!-------------------------------------------------------------------
47      !!                ***  ROUTINE lim_thd_sal  ***       
48      !! ** Purpose :
49      !!        This routine computes new salinities in the ice
50      !!
51      !! ** Method  :  4 possibilities
52      !!               -> num_sal = 1 -> constant salinity for z,t
53      !!               -> num_sal = 2 -> S = S(z,t) [simple Vancoppenolle et al 2005]
54      !!               -> num_sal = 3 -> S = S(z)   [multiyear ice]
55      !!               -> num_sal = 4 -> S = S(h)   [Cox and Weeks 74]
56      !!           
57      !! ** Steps
58      !!
59      !! ** Arguments
60      !!
61      !! ** Inputs / Outputs
62      !!
63      !! ** External
64      !!
65      !! ** References
66      !!
67      !! ** History  :
68      !!
69      !! "Je ne suis pour l'instant qu'a 80% de ma condition, mais c'est
70      !!  les 30% qui restent qui seront les plus difficiles"
71      !!                                           E. Mpenza
72      !!
73      !!-------------------------------------------------------------------
74      !! History :
75      !!   ori  !  03-05 M. Vancoppenolle UCL-ASTR first coding for LIM-1D
76      !!   3.0  !  05-12 Routine rewritten for the 3-D version
77      !!---------------------------------------------------------------------
78      !!
79      !! * Local variables
80      INTEGER, INTENT(in) :: &
81         kideb, kiut             !: thickness category index
82
83      INTEGER ::             &
84         ji, jk ,            &   !: geographic and layer index
85         zji, zjj
86
87      REAL(wp) ::            &
88         zsold,              &   !: old salinity
89         zeps=1.0e-06   ,    &   !: very small
90         iflush         ,    &   !: flushing (1) or not (0)
91         iaccrbo        ,    &   !: bottom accretion (1) or not (0)
92         igravdr        ,    &   !: gravity drainage or not
93         isnowic        ,    &   !: gravity drainage or not
94         i_ice_switch   ,    &   !: ice thickness above a certain treshold or not
95         ztmelts        ,    &   !: freezing point of sea ice
96         zaaa           ,    &   !: dummy factor
97         zbbb           ,    &   !: dummy factor
98         zccc           ,    &   !: dummy factor
99         zdiscrim                !: dummy factor
100
101      REAL(wp), DIMENSION(jpij) ::          &
102         ze_init        ,    &   !initial total enthalpy
103         zhiold         ,    &
104         zsiold
105
106      !!---------------------------------------------------------------------
107
108      !------------------------------------------------------------------------------|
109      ! 1) Constant salinity, constant in time                                       |
110      !------------------------------------------------------------------------------|
111
112      IF (num_sal.eq.1) THEN
113
114         !         WRITE(numout,*)
115         !         WRITE(numout,*) 'lim_thd_sal : Ice salinity computation module ', &
116         !         num_sal
117         !         WRITE(numout,*) '~~~~~~~~~~~~'
118
119         DO jk = 1, nlay_i
120            DO ji = kideb, kiut
121               s_i_b(ji,jk) =  bulk_sal
122            END DO ! ji
123         END DO ! jk
124
125         DO ji = kideb, kiut
126            sm_i_b(ji)      =  bulk_sal 
127         END DO ! ji
128
129      ENDIF ! num_sal .EQ. 1
130
131      !------------------------------------------------------------------------------|
132      !  Module 2 : Constant salinity varying in time                                |
133      !------------------------------------------------------------------------------|
134
135      IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) THEN
136
137         !         WRITE(numout,*)
138         !         WRITE(numout,*) 'lim_thd_sal : Ice salinity computation module ', &
139         !         num_sal
140         !         WRITE(numout,*) '~~~~~~~~~~~'
141         !         WRITE(numout,*)
142
143         !---------------------------------
144         ! Thickness at previous time step
145         !---------------------------------
146         DO ji = kideb, kiut
147            zhiold(ji)   =  ht_i_b(ji) - dh_i_bott(ji) - dh_snowice(ji) -     &
148               dh_i_surf(ji)
149         END DO ! ji
150
151         !---------------------
152         ! Global heat content
153         !---------------------
154
155         ze_init(:)  =  0.0
156         DO jk = 1, nlay_i
157            DO ji = kideb, kiut
158               ze_init(ji) = ze_init(ji) + q_i_b(ji,jk) * ht_i_b(ji) / nlay_i
159            END DO ! ji
160         END DO ! jk
161
162         DO ji = kideb, kiut
163
164            !----------
165            ! Switches
166            !----------
167
168            ! iflush  : 1 if summer
169            iflush       =  MAX( 0.0 , SIGN ( 1.0 , t_su_b(ji) - rtt ) ) 
170            ! igravdr : 1 if t_su lt t_bo
171            igravdr      =  MAX( 0.0 , SIGN ( 1.0 , t_bo_b(ji) - t_su_b(ji) ) )
172            ! iaccrbo : 1 if bottom accretion
173            iaccrbo      =  MAX( 0.0 , SIGN ( 1.0 , dh_i_bott(ji) ) )
174            ! isnowic : 1 if snow ice formation
175            i_ice_switch = 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_i_b(ji) + 1.0e-2 ) )
176            isnowic      = 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - dh_snowice(ji) ) ) * &
177               i_ice_switch
178
179            !---------------------
180            ! Salinity tendencies
181            !---------------------
182
183            ! drainage by gravity drainage
184            dsm_i_gd_1d(ji) = - igravdr *                                     & 
185               MAX( sm_i_b(ji) - sal_G , 0.0 ) /             &
186               time_G * rdt_ice 
187
188            ! drainage by flushing 
189            dsm_i_fl_1d(ji)  = - iflush *                                     &
190               MAX( sm_i_b(ji) - sal_F , 0.0 ) /             & 
191               time_F * rdt_ice
192
193            !-----------------
194            ! Update salinity   
195            !-----------------
196
197            ! only drainage terms ( gravity drainage and flushing )
198            ! snow ice / bottom sources are added in lim_thd_ent
199            ! to conserve energy
200            zsiold(ji) = sm_i_b(ji)
201            sm_i_b(ji) = sm_i_b(ji)                                           &
202               + dsm_i_fl_1d(ji) + dsm_i_gd_1d(ji)  !                 &
203
204            ! if no ice, salinity eq 0.1
205            i_ice_switch  = 1.0 - MAX ( 0.0, SIGN (1.0 , - ht_i_b(ji) ) )
206            sm_i_b(ji)    = i_ice_switch*sm_i_b(ji) + s_i_min * ( 1.0 -       &
207               i_ice_switch )
208         END DO ! ji
209
210         ! Salinity profile
211         CALL lim_var_salprof1d(kideb,kiut)
212
213         !----------------------------
214         ! Heat flux - brine drainage
215         !----------------------------
216
217         DO ji = kideb, kiut
218            ! iflush  : 1 if summer
219            iflush       =  MAX( 0.0 , SIGN ( 1.0 , t_su_b(ji) - rtt ) ) 
220            ! igravdr : 1 if t_su lt t_bo
221            igravdr      =  MAX( 0.0 , SIGN ( 1.0 , t_bo_b(ji) - t_su_b(ji) ) ) 
222            ! iaccrbo : 1 if bottom accretion
223            iaccrbo      =  MAX( 0.0 , SIGN ( 1.0 , dh_i_bott(ji) ) )
224
225            fhbri_1d(ji) = 0.0
226         END DO ! ji
227
228         !----------------------------
229         ! Salt flux - brine drainage
230         !----------------------------
231         DO ji = kideb, kiut
232            i_ice_switch  = 1.0 - MAX ( 0.0, SIGN (1.0 , - ht_i_b(ji) ) )
233            fsbri_1d(ji) = fsbri_1d(ji) -  &
234               i_ice_switch * rhoic * a_i_b(ji) * ht_i_b(ji) *  &
235               ( MAX(dsm_i_gd_1d(ji) + dsm_i_fl_1d(ji), &
236               sm_i_b(ji) - zsiold(ji) ) ) / rdt_ice
237            IF ( num_sal .EQ. 4 ) fsbri_1d(ji) = 0.0
238
239         END DO ! ji
240
241         ! Only necessary for conservation check since salinity is modified
242         !--------------------
243         ! Temperature update
244         !--------------------
245         DO jk = 1, nlay_i
246
247            DO ji = kideb, kiut
248
249               ztmelts    =  -tmut*s_i_b(ji,jk) + rtt
250               !Conversion q(S,T) -> T (second order equation)
251               zaaa         =  cpic
252               zbbb         =  ( rcp - cpic ) * ( ztmelts - rtt ) + &
253                  q_i_b(ji,jk) / rhoic - lfus
254               zccc         =  lfus * ( ztmelts - rtt )
255               zdiscrim     =  SQRT( MAX(zbbb*zbbb - 4.0*zaaa*zccc,0.0) )
256               t_i_b(ji,jk) =  rtt - ( zbbb + zdiscrim ) / & 
257                  ( 2.0 *zaaa )
258            END DO !ji
259
260         END DO !jk
261
262      ENDIF ! num_sal .EQ. 2
263
264      !------------------------------------------------------------------------------|
265      !  Module 3 : Profile of salinity, constant in time                            |
266      !------------------------------------------------------------------------------|
267
268      IF ( num_sal .EQ. 3 ) THEN
269
270         WRITE(numout,*)
271         WRITE(numout,*) 'lim_thd_sal : Ice salinity computation module ', &
272            num_sal
273         WRITE(numout,*) '~~~~~~~~~~~~'
274
275         CALL lim_var_salprof1d(kideb,kiut)
276
277      ENDIF ! num_sal .EQ. 3
278
279      !------------------------------------------------------------------------------|
280      !  Module 4 : Constant salinity varying in time                                |
281      !------------------------------------------------------------------------------|
282
283      ! Cox and Weeks, 1974
284      IF (num_sal.eq.5) THEN
285
286         WRITE(numout,*)
287         WRITE(numout,*) 'lim_thd_sal : Ice salinity computation module ', &
288            num_sal
289         WRITE(numout,*) '~~~~~~~~~~~~'
290
291         DO ji = kideb, kiut
292
293            zsold = sm_i_b(ji)
294
295            IF (ht_i_b(ji).lt.0.4) THEN
296               sm_i_b(ji)    = 14.24 - 19.39*ht_i_b(ji) 
297            ELSE
298               sm_i_b(ji)    =  7.88 - 1.59*ht_i_b(ji)
299               sm_i_b(ji)    = MIN(sm_i_b(ji),zsold) 
300            ENDIF
301
302            IF ( ht_i_b(ji) .GT. 3.06918239 ) THEN
303               sm_i_b(ji)     = 3.0
304            ENDIF
305
306            DO jk = 1, nlay_i
307               s_i_b(ji,jk)   = sm_i_b(ji)
308            END DO
309
310         END DO ! ji
311
312      ENDIF ! num_sal
313
314      !------------------------------------------------------------------------------|
315      ! 5) Computation of salt flux due to Bottom growth
316      !------------------------------------------------------------------------------|
317
318      IF ( num_sal .EQ. 4 ) THEN
319         DO ji = kideb, kiut
320            zji                 = MOD( npb(ji) - 1, jpi ) + 1
321            zjj                 = ( npb(ji) - 1 ) / jpi + 1
322            fseqv_1d(ji) = fseqv_1d(ji)              + & 
323               ( sss_m(zji,zjj) - bulk_sal    ) * & 
324               rhoic * a_i_b(ji) * &
325               MAX( dh_i_bott(ji) , 0.0 ) / rdt_ice
326         END DO
327      ELSE
328         DO ji = kideb, kiut
329            zji                 = MOD( npb(ji) - 1, jpi ) + 1
330            zjj                 = ( npb(ji) - 1 ) / jpi + 1
331            fseqv_1d(ji) = fseqv_1d(ji)              + & 
332               ( sss_m(zji,zjj) - s_i_new(ji) ) * & 
333               rhoic * a_i_b(ji) * &
334               MAX( dh_i_bott(ji) , 0.0 ) / rdt_ice
335         END DO ! ji
336      ENDIF
337
338      !-- End of salinity computations
339   END SUBROUTINE lim_thd_sal
340   !==============================================================================
341
342   SUBROUTINE lim_thd_sal_init
343      !!-------------------------------------------------------------------
344      !!                  ***  ROUTINE lim_thd_sal_init  ***
345      !!
346      !! ** Purpose :   initialization of ice salinity parameters
347      !!
348      !! ** Method  : Read the namicesal namelist and check the parameter
349      !!       values called at the first timestep (nit000)
350      !!
351      !! ** input   :   Namelist namicesal
352      !!
353      !! history :
354      !!   3.0  !  July 2005 M. Vancoppenolle  Original code
355      !!-------------------------------------------------------------------
356      NAMELIST/namicesal/ num_sal, bulk_sal, sal_G, time_G, sal_F, time_F, &
357         s_i_max, s_i_min, s_i_0, s_i_1
358      !!-------------------------------------------------------------------
359
360      ! Read Namelist namicesal
361      REWIND ( numnam_ice )
362      READ   ( numnam_ice  , namicesal )
363      IF(lwp) THEN
364         WRITE(numout,*)
365         WRITE(numout,*) 'lim_thd_sal_init : Ice parameters for salinity '
366         WRITE(numout,*) '~~~~~~~~~~~~~~~~'
367         WRITE(numout,*) ' switch for salinity num_sal        : ', num_sal
368         WRITE(numout,*) ' bulk salinity value if num_sal = 1 : ', bulk_sal
369         WRITE(numout,*) ' restoring salinity for GD          : ', sal_G
370         WRITE(numout,*) ' restoring time for GD              : ', time_G
371         WRITE(numout,*) ' restoring salinity for flushing    : ', sal_F
372         WRITE(numout,*) ' restoring time for flushing        : ', time_F
373         WRITE(numout,*) ' Maximum tolerated ice salinity     : ', s_i_max
374         WRITE(numout,*) ' Minimum tolerated ice salinity     : ', s_i_min
375         WRITE(numout,*) ' 1st salinity for salinity profile  : ', s_i_0
376         WRITE(numout,*) ' 2nd salinity for salinity profile  : ', s_i_1
377      ENDIF
378
379   END SUBROUTINE lim_thd_sal_init
380
381#else
382   !!----------------------------------------------------------------------
383   !!   Default option         Empty Module                No sea-ice model
384   !!----------------------------------------------------------------------
385CONTAINS
386   SUBROUTINE lim_thd_sal        ! Empty routine
387   END SUBROUTINE lim_thd_sal
388#endif
389   !!======================================================================
390END MODULE limthd_sal
Note: See TracBrowser for help on using the repository browser.