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

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90 @ 2287

Last change on this file since 2287 was 2287, checked in by smasson, 14 years ago

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