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 @ 1146

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

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