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

source: trunk/NEMO/LIM_SRC_3/limthd_sal.F90 @ 1156

Last change on this file since 1156 was 1156, checked in by rblod, 16 years ago

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