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

source: branches/dev_002_LIM/NEMO/LIM_SRC_3/limthd_sal.F90 @ 825

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

dev_002_LIM : add the LIM 3.0 component, see ticketr: #71
File size: 18.8 KB
Line 
1MODULE limthd_sal
2#if defined key_lim3
3   !!======================================================================
4   !!                       ***  MODULE limthd_sal ***
5   !!                   computation salinity variations in
6   !!                               the ice
7   !!======================================================================
8
9   !!----------------------------------------------------------------------
10   !!   lim_thd_sal : salinity variations in the ice
11   !! * Modules used
12   USE par_oce          ! ocean parameters
13   USE phycst           ! physical constants (ocean directory)
14   USE ice_oce          ! ice variables
15   USE thd_ice
16   USE iceini
17   USE ice
18   USE limistate
19   USE in_out_manager
20   USE limvar
21   USE par_ice
22   USE limicepoints
23
24   IMPLICIT NONE
25   PRIVATE
26
27   !! * Routine accessibility
28   PUBLIC lim_thd_sal        ! 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   CONTAINS
38
39   SUBROUTINE lim_thd_sal(kideb,kiut,jl)
40      !!-------------------------------------------------------------------
41      !!                ***  ROUTINE lim_thd_sal  ***       
42      !! ** Purpose :
43      !!        This routine computes new salinities in the ice
44      !!
45      !! ** Method  :  4 possibilities
46      !!               -> num_sal = 1 -> constant salinity for z,t
47      !!               -> num_sal = 2 -> S = S(z,t) [simple Vancoppenolle et al 2005]
48      !!               -> num_sal = 3 -> S = S(z)   [multiyear ice]
49      !!               -> num_sal = 4 -> S = S(h)   [Cox and Weeks 74]
50      !!           
51      !! ** Steps
52      !!
53      !! ** Arguments
54      !!
55      !! ** Inputs / Outputs
56      !!
57      !! ** External
58      !!
59      !! ** References
60      !!
61      !! ** History  :
62      !!-------------------------------------------------------------------
63      !! History :
64      !!   ori  !  03-05 M. Vancoppenolle UCL-ASTR first coding for LIM-1D
65      !!   3.0  !  05-12 Routine rewritten for the 3-D version
66      !!---------------------------------------------------------------------
67      !! * Local variables
68      INTEGER, INTENT(in) :: &
69         kideb, kiut, jl         !: thickness category index
70
71      INTEGER ::             &
72         ji, jk ,            &   !: geographic and layer index
73         zji, zjj
74
75      REAL(wp) ::            &
76!        ds_i_gd,          &     !: salt loss by gravity drainage
77!        ds_i_seg,         &     !: salt gain by segregation of seawater at the bottom
78!        ds_i_flush,       &     !: summer salt loss by flushing
79         zargtemp,           &   !: used to compute Scharzacher profile
80         zsold,              &   !: old salinity
81!        ds_i_snice,       &     !: increase of salinity due to snow ice
82         zeps=1.0e-06   ,    &   !: very small
83         iflush         ,    &   !: flushing (1) or not (0)
84         iaccrbo        ,    &   !: bottom accretion (1) or not (0)
85         igravdr        ,    &   !: gravity drainage or not
86         isnowic        ,    &   !: gravity drainage or not
87         i_ice_switch   ,    &   !: ice thickness above a certain treshold or not
88         ztmelts        ,    &   !: freezing point of sea ice
89         zs_snowice     ,    &   !: salinity of forming snow ice
90         zaaa           ,    &   !: dummy factor
91         zbbb           ,    &   !: dummy factor
92         zccc           ,    &   !: dummy factor
93         zdiscrim                !: dummy factor
94 
95      REAL(wp), DIMENSION(jpij) ::          &
96         ze_init        ,    &   !initial total enthalpy
97         ze_end         ,    &   !final total enthalpy
98         zhiold         ,    &
99         zsiold
100
101      !!---------------------------------------------------------------------
102
103      IF ( ( numit == nstart ) .AND. ( jl == 1 ) ) &
104         CALL lim_thd_sal_init   ! Initialization
105
106!------------------------------------------------------------------------------|
107! 1) Constant salinity, constant in time                                       |
108!------------------------------------------------------------------------------|
109
110      IF (num_sal.eq.1) THEN
111
112         WRITE(numout,*)
113         WRITE(numout,*) 'lim_thd_sal : Ice salinity computation module ', &
114         num_sal
115         WRITE(numout,*) '~~~~~~~~~~~~'
116
117         DO jk = 1, nlay_i
118            DO ji = kideb, kiut
119               s_i_b(ji,jk) =  bulk_sal
120            END DO ! ji
121         END DO ! jk
122
123         DO ji = kideb, kiut
124            sm_i_b(ji)      =  bulk_sal 
125         END DO ! ji
126
127      ENDIF ! num_sal .EQ. 1
128
129!------------------------------------------------------------------------------|
130!  Module 2 : Constant salinity varying in time                                |
131!------------------------------------------------------------------------------|
132
133      IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) THEN
134
135         WRITE(numout,*)
136         WRITE(numout,*) 'lim_thd_sal : Ice salinity computation module ', &
137         num_sal
138         WRITE(numout,*) '~~~~~~~~~~~'
139         WRITE(numout,*)
140
141         !---------------------------------
142         ! Thickness at previous time step
143         !---------------------------------
144         DO ji = kideb, kiut
145            zhiold(ji)   =  ht_i_b(ji) - dh_i_bott(ji) - dh_snowice(ji) -     &
146                            dh_i_surf(ji)
147         END DO ! ji
148
149         !---------------------
150         ! Global heat content
151         !---------------------
152
153         ze_init(:)  =  0.0
154         DO jk = 1, nlay_i
155            DO ji = kideb, kiut
156               ze_init(ji) = ze_init(ji) + q_i_b(ji,jk) * ht_i_b(ji) / nlay_i
157            END DO ! ji
158         END DO ! jk
159
160 !       IF (jiindex_1d .GT. 0 ) THEN
161 !       WRITE(numout,*) ' s_i_b   : ', s_i_b(jiindex_1d,1:nlay_i)
162 !       WRITE(numout,*) ' t_i_b   : ', t_i_b(jiindex_1d,1:nlay_i)
163 !       WRITE(numout,*) ' q_i_b   : ', q_i_b(jiindex_1d,1:nlay_i)
164 !       ENDIF
165
166         DO ji = kideb, kiut
167
168            !----------
169            ! Switches
170            !----------
171
172            ! iflush  : 1 if summer
173            iflush       =  MAX( 0.0 , SIGN ( 1.0 , t_su_b(ji) - rtt ) ) 
174            ! igravdr : 1 if t_su lt t_bo
175            igravdr      =  MAX( 0.0 , SIGN ( 1.0 , t_bo_b(ji) - t_su_b(ji) ) )
176            ! iaccrbo : 1 if bottom accretion
177            iaccrbo      =  MAX( 0.0 , SIGN ( 1.0 , dh_i_bott(ji) ) )
178            ! isnowic : 1 if snow ice formation
179            i_ice_switch = 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_i_b(ji) + 1.0d-2 ) )
180            isnowic      = 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - dh_snowice(ji) ) ) * &
181                           i_ice_switch
182
183            !---------------------
184            ! Salinity tendencies
185            !---------------------
186
187            ! drainage by gravity drainage
188            dsm_i_gd_1d(ji) = - igravdr *                                     & 
189                                MAX( sm_i_b(ji) - sal_G , 0.0 ) /             &
190                                time_G * rdt_ice 
191
192            ! drainage by flushing 
193            dsm_i_fl_1d(ji)  = - iflush *                                     &
194                                MAX( sm_i_b(ji) - sal_F , 0.0 ) /             & 
195                                time_F * rdt_ice
196
197            !-----------------
198            ! Update salinity   
199            !-----------------
200
201            ! only drainage terms ( gravity drainage and flushing )
202            ! snow ice / bottom sources are added in lim_thd_ent
203            ! to conserve energy
204            zsiold(ji) = sm_i_b(ji)
205            sm_i_b(ji) = sm_i_b(ji)                                           &
206                       + dsm_i_fl_1d(ji) + dsm_i_gd_1d(ji)  !                 &
207
208            ! if no ice, salinity eq 0.1
209            i_ice_switch  = 1.0 - MAX ( 0.0, SIGN (1.0 , - ht_i_b(ji) ) )
210            sm_i_b(ji)    = i_ice_switch*sm_i_b(ji) + s_i_min * ( 1.0 -       &
211                            i_ice_switch )
212         END DO ! ji
213
214         ! Salinity profile
215         CALL lim_var_salprof1d(kideb,kiut)
216
217         !--------------------------------------
218         ! Energy of melting for each ice layer
219         !--------------------------------------
220         ! q=q(S,T) J.m-3
221!        DO jk = 1, nlay_i
222!           DO ji = kideb, kiut
223!           ! New heat content after desalination
224!              ztmelts      = - tmut * s_i_b(ji,jk) + rtt ! Tfi in K
225!              q_i_b(ji,jk) = rhoic *                                         &
226!                           ( cpic * ( ztmelts-t_i_b(ji,jk) )                 &
227!                           + lfus * ( 1.0 - (ztmelts-rtt) /                  &
228!                             MIN( ( t_i_b(ji,jk) - rtt ) , -zeps ) )         &
229!                           - rcp  * ( ztmelts-rtt ) )
230!           END DO ! ji
231!        END DO
232
233         !--------------------
234         ! Total heat content
235         !--------------------
236!        ze_end(:) = 0.0
237!        DO jk = 1, nlay_i
238!           DO ji = kideb, kiut
239!              ze_end(ji)  = ze_end(ji) + q_i_b(ji,jk) * ht_i_b(ji) / nlay_i
240!           END DO
241!        END DO
242           
243!        IF (jiindex_1d .GT. 0 ) THEN
244!           WRITE(numout,*) ' ze_init : ', ze_init(jiindex_1d) / rdt_ice
245!           WRITE(numout,*) ' ze_end  : ', ze_end (jiindex_1d) /  rdt_ice
246!           WRITE(numout,*) ' s_i_b   : ', s_i_b(jiindex_1d,1:nlay_i)
247!           WRITE(numout,*) ' t_i_b   : ', t_i_b(jiindex_1d,1:nlay_i)
248!           WRITE(numout,*) ' q_i_b   : ', q_i_b(jiindex_1d,1:nlay_i)
249!        ENDIF
250
251         !----------------------------
252         ! Heat flux - brine drainage
253         !----------------------------
254
255!        WRITE(numout,*) ' jiindex, jjindex : ', jiindex, jjindex
256!        WRITE(numout,*) ' jiindex_1d : ', jiindex_1d
257!        WRITE(numout,*) ' kideb, kiut: ', kideb, kiut
258
259         DO ji = kideb, kiut
260            ! iflush  : 1 if summer
261            iflush       =  MAX( 0.0 , SIGN ( 1.0 , t_su_b(ji) - rtt ) ) 
262            ! igravdr : 1 if t_su lt t_bo
263            igravdr      =  MAX( 0.0 , SIGN ( 1.0 , t_bo_b(ji) - t_su_b(ji) ) ) 
264            ! iaccrbo : 1 if bottom accretion
265            iaccrbo      =  MAX( 0.0 , SIGN ( 1.0 , dh_i_bott(ji) ) )
266
267         !  IF ( ( iflush .EQ. 1 )  .OR. ( iaccrbo .EQ. 1 ) ) THEN
268!           fhbri_1d(ji) = fhbri_1d(ji) + ( ze_end(ji) - ze_init(ji) ) * &
269!                          a_i_b(ji) / at_i_b(ji) / rdt_ice
270                           ! ze_end GT ze_init... ?
271         !  ENDIF
272            fhbri_1d(ji) = 0.0
273
274            !++++++
275!           zji                 = MOD( npb(ji) - 1, jpi ) + 1
276!           zjj                 = ( npb(ji) - 1 ) / jpi + 1
277!           fhbricat(zji,zjj,jl)  =   ( ze_end(ji) - ze_init(ji) ) / &
278!                                   rdt_ice
279            !++++++
280!           ftotal_fin(ji) = ftotal_fin(ji) + fhbricat(zji,zjj,jl)
281!           IF ( ji .EQ. jiindex_1d) WRITE(numout,*) ' fhbri : ', &
282!                                  fhbricat(zji,zjj,jl) * rdt_ice
283!           IF ( ji .EQ. jiindex_1d ) &
284!              WRITE(numout,*) ' ftotal_fin : ', ftotal_fin(jiindex_1d)
285!           !+++++
286
287         END DO ! ji
288
289!        IF (jiindex_1d .GT. 0 ) THEN
290!           WRITE(numout,*) ' Category : ', jl
291!           WRITE(numout,*) ' jiindex_1d : ', jiindex_1d
292!           WRITE(numout,*) ' fhbricat : ', fhbricat(jiindex,jjindex,jl)
293!           WRITE(numout,*) ' fhbri_1d : ', fhbri_1d(jiindex_1d)
294!           WRITE(numout,*) ' ze_end   : ', ze_end(jiindex_1d)
295!           WRITE(numout,*) ' ze_init  : ', ze_init(jiindex_1d)
296!           WRITE(numout,*) ' a_i_b    : ', a_i_b(jiindex_1d)
297!           WRITE(numout,*) ' at_i_b   : ', at_i_b(jiindex_1d)
298!           WRITE(numout,*) ' rdt_ice  : ', rdt_ice
299!        ENDIF
300
301         !----------------------------
302         ! Salt flux - brine drainage
303         !----------------------------
304         DO ji = kideb, kiut
305            i_ice_switch  = 1.0 - MAX ( 0.0, SIGN (1.0 , - ht_i_b(ji) ) )
306            fsbri_1d(ji) = fsbri_1d(ji) -  &
307                           i_ice_switch * rhoic * a_i_b(ji) * ht_i_b(ji) *  &
308                           ( MAX(dsm_i_gd_1d(ji) + dsm_i_fl_1d(ji), &
309                             sm_i_b(ji) - zsiold(ji) ) ) / rdt_ice
310            IF ( num_sal .EQ. 4 ) fsbri_1d(ji) = 0.0
311
312         END DO ! ji
313
314         ! Only necessary for conservation check since salinity is modified
315         !--------------------
316         ! Temperature update
317         !--------------------
318         DO jk = 1, nlay_i
319
320            DO ji = kideb, kiut
321
322               ztmelts    =  -tmut*s_i_b(ji,jk) + rtt
323               !Conversion q(S,T) -> T (second order equation)
324               zaaa         =  cpic
325               zbbb         =  ( rcp - cpic ) * ( ztmelts - rtt ) + &
326                               q_i_b(ji,jk) / rhoic - lfus
327               zccc         =  lfus * ( ztmelts - rtt )
328               zdiscrim     =  SQRT( MAX(zbbb*zbbb - 4.0*zaaa*zccc,0.0) )
329               t_i_b(ji,jk) =  rtt - ( zbbb + zdiscrim ) / & 
330                                     ( 2.0 *zaaa )
331            END DO !ji
332
333         END DO !jk
334
335      ENDIF ! num_sal .EQ. 2
336
337!------------------------------------------------------------------------------|
338!  Module 3 : Profile of salinity, constant in time                            |
339!------------------------------------------------------------------------------|
340
341      IF ( num_sal .EQ. 3 ) THEN
342
343         WRITE(numout,*)
344         WRITE(numout,*) 'lim_thd_sal : Ice salinity computation module ', &
345         num_sal
346         WRITE(numout,*) '~~~~~~~~~~~~'
347
348         CALL lim_var_salprof1d(kideb,kiut)
349
350      ENDIF ! num_sal .EQ. 3
351
352!------------------------------------------------------------------------------|
353!  Module 4 : Constant salinity varying in time                                |
354!------------------------------------------------------------------------------|
355
356      ! Cox and Weeks, 1974
357      IF (num_sal.eq.5) THEN
358
359         WRITE(numout,*)
360         WRITE(numout,*) 'lim_thd_sal : Ice salinity computation module ', &
361         num_sal
362         WRITE(numout,*) '~~~~~~~~~~~~'
363
364         DO ji = kideb, kiut
365
366            zsold = sm_i_b(ji)
367
368            IF (ht_i_b(ji).lt.0.4) THEN
369               sm_i_b(ji)    = 14.24 - 19.39*ht_i_b(ji) 
370            ELSE
371               sm_i_b(ji)    =  7.88 - 1.59*ht_i_b(ji)
372               sm_i_b(ji)    = MIN(sm_i_b(ji),zsold) 
373            ENDIF
374         
375            IF ( ht_i_b(ji) .GT. 3.06918239 ) THEN
376               sm_i_b(ji)     = 3.0
377            ENDIF
378
379      !!--multiyear sea-ice salinity
380      !!-- ca va pas alors on vire bordel de djeuille
381!            if (ht.gt.360) then
382!           do jk = 1, nlay_i
383!           s_i_b(ji,jk) = 1.58 + 0.18*ht_i_b(ji)
384!           end do
385!        endif
386            DO jk = 1, nlay_i
387               s_i_b(ji,jk)   = sm_i_b(ji)
388            END DO
389     
390         END DO ! ji
391
392      ENDIF ! num_sal
393
394!------------------------------------------------------------------------------|
395! 5) Computation of salt flux due to Bottom growth
396!------------------------------------------------------------------------------|
397
398      IF ( num_sal .EQ. 4 ) THEN
399         DO ji = kideb, kiut
400            zji                 = MOD( npb(ji) - 1, jpi ) + 1
401            zjj                 = ( npb(ji) - 1 ) / jpi + 1
402            fseqv_1d(ji) = fseqv_1d(ji)              + & 
403                           ( sss_io(zji,zjj) - bulk_sal    ) * & 
404                           rhoic * a_i_b(ji) * &
405                           MAX( dh_i_bott(ji) , 0.0 ) / rdt_ice
406         END DO
407      ELSE
408         DO ji = kideb, kiut
409            zji                 = MOD( npb(ji) - 1, jpi ) + 1
410            zjj                 = ( npb(ji) - 1 ) / jpi + 1
411            fseqv_1d(ji) = fseqv_1d(ji)              + & 
412                           ( sss_io(zji,zjj) - s_i_new(ji) ) * & 
413                             rhoic * a_i_b(ji) * &
414                             MAX( dh_i_bott(ji) , 0.0 ) / rdt_ice
415         END DO ! ji
416      ENDIF
417
418!+++++
419!           IF ( (zji.EQ.jiindex) .AND. (zjj.EQ.jjindex) ) THEN
420!              WRITE(numout,*) ' *** FS_EQV *** '
421!              WRITE(numout,*) ' limthd_sal '
422!              WRITE(numout,*) ' fseqv      ', fseqv_1d(ji)
423!              WRITE(numout,*) ' sss_io     ', sss_io(zji,zjj)
424!              WRITE(numout,*) ' s_i_new    ', s_i_new(ji)
425!              WRITE(numout,*) ' dh_i_bott  ', dh_i_bott(ji)
426!              WRITE(numout,*) ' a_i_b      ', a_i_b(ji)
427!              WRITE(numout,*) ' *** FS_BRI *** '
428!              WRITE(numout,*) ' fsbri      ', fsbri_1d(ji)
429!              WRITE(numout,*) ' ht_i       ', ht_i_b(ji)
430!              WRITE(numout,*) ' dsm_i_gd   ', dsm_i_gd_1d(ji)
431!              WRITE(numout,*) ' dsm_i_fl   ', dsm_i_fl_1d(ji)
432!           ENDIF
433
434!+++++
435
436!-- End of salinity computations
437   END SUBROUTINE lim_thd_sal
438!==============================================================================
439
440   SUBROUTINE lim_thd_sal_init
441      !!-------------------------------------------------------------------
442      !!                  ***  ROUTINE lim_thd_sal_init  ***
443      !!
444      !! ** Purpose :   initialization of ice salinity parameters
445      !!
446      !! ** Method  : Read the namicesal namelist and check the parameter
447      !!       values called at the first timestep (nit000)
448      !!
449      !! ** input   :   Namelist namicesal
450      !!
451      !! history :
452      !!   3.0  !  July 2005 M. Vancoppenolle  Original code
453      !!-------------------------------------------------------------------
454      NAMELIST/namicesal/ num_sal, bulk_sal, sal_G, time_G, sal_F, time_F, &
455                          s_i_max, s_i_min, s_i_0, s_i_1
456      !!-------------------------------------------------------------------
457
458      ! Read Namelist namicesal
459      REWIND ( numnam_ice )
460      READ   ( numnam_ice  , namicesal )
461      IF(lwp) THEN
462         WRITE(numout,*)
463         WRITE(numout,*) 'lim_thd_sal_init : Ice parameters for salinity '
464         WRITE(numout,*) '~~~~~~~~~~~~~~~~'
465         WRITE(numout,*) ' switch for salinity num_sal        : ', num_sal
466         WRITE(numout,*) ' bulk salinity value if num_sal = 1 : ', bulk_sal
467         WRITE(numout,*) ' restoring salinity for GD          : ', sal_G
468         WRITE(numout,*) ' restoring time for GD              : ', time_G
469         WRITE(numout,*) ' restoring salinity for flushing    : ', sal_F
470         WRITE(numout,*) ' restoring time for flushing        : ', time_F
471         WRITE(numout,*) ' Maximum tolerated ice salinity     : ', s_i_max
472         WRITE(numout,*) ' Minimum tolerated ice salinity     : ', s_i_min
473         WRITE(numout,*) ' 1st salinity for salinity profile  : ', s_i_0
474         WRITE(numout,*) ' 2nd salinity for salinity profile  : ', s_i_1
475      ENDIF
476
477   END SUBROUTINE lim_thd_sal_init
478
479#else
480   !!----------------------------------------------------------------------
481   !!   Default option         Empty Module                No sea-ice model
482   !!----------------------------------------------------------------------
483CONTAINS
484   SUBROUTINE lim_thd_sal        ! Empty routine
485   END SUBROUTINE lim_thd_sal
486#endif
487   !!======================================================================
488END MODULE limthd_sal
Note: See TracBrowser for help on using the repository browser.