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

Last change on this file since 867 was 842, checked in by ctlod, 16 years ago

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