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

source: branches/2011/dev_r2787_LOCEAN3_TRA_TRP/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90 @ 2833

Last change on this file since 2833 was 2777, checked in by smasson, 13 years ago

LIM3 compiling and (partly?) running in v3_3_1, see ticket#817
  • Property svn:keywords set to Id
File size: 12.8 KB
Line 
1MODULE limthd_sal
2   !!======================================================================
3   !!                       ***  MODULE limthd_sal ***
4   !! LIM-3 sea-ice :  computation of salinity variations in the ice
5   !!======================================================================
6   !! History :   -   ! 2003-05 (M. Vancoppenolle) UCL-ASTR first coding for LIM3-1D
7   !!            3.0  ! 2005-12 (M. Vancoppenolle) adapted to the 3-D version
8   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation
9   !!---------------------------------------------------------------------
10#if defined key_lim3
11   !!----------------------------------------------------------------------
12   !!   'key_lim3'                                      LIM-3 sea-ice model
13   !!----------------------------------------------------------------------
14   !!   lim_thd_sal : salinity variations in the ice
15   !!----------------------------------------------------------------------
16   USE par_oce          ! ocean parameters
17   USE phycst           ! physical constants (ocean directory)
18   USE sbc_oce          ! Surface boundary condition: ocean fields
19   USE ice              ! LIM variables
20   USE par_ice          ! LIM parameters
21   USE thd_ice          ! LIM thermodynamics
22   USE limvar           ! LIM variables
23   USE in_out_manager   ! I/O manager
24   USE lib_mpp         ! MPP library
25
26   IMPLICIT NONE
27   PRIVATE
28
29   PUBLIC   lim_thd_sal        ! called by limthd module
30   PUBLIC   lim_thd_sal_init   ! called by iceini module
31
32   !!----------------------------------------------------------------------
33   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
34   !! $Id$
35   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
36   !!----------------------------------------------------------------------
37CONTAINS
38
39   SUBROUTINE lim_thd_sal( kideb, kiut )
40      !!-------------------------------------------------------------------
41      !!                ***  ROUTINE lim_thd_sal  ***   
42      !!   
43      !! ** Purpose :   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      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
52      USE wrk_nemo, ONLY:   wrk_1d_1, wrk_1d_2, wrk_1d_3
53      !
54      INTEGER, INTENT(in) ::  kideb, kiut   ! thickness category index
55      !
56      INTEGER  ::   ji, jk     ! dummy loop indices
57      INTEGER  ::   zji, zjj   ! local integers
58      REAL(wp) ::   zsold, iflush, iaccrbo, igravdr, isnowic, i_ice_switch,  ztmelts   ! local scalars
59      REAL(wp) ::   zaaa, zbbb, zccc, zdiscrim   ! local scalars
60      !
61      REAL(wp), POINTER, DIMENSION(:) ::   ze_init, zhiold, zsiold
62      !!---------------------------------------------------------------------
63
64      IF(  wrk_in_use(1, 1,2,3)  ) THEN
65         CALL ctl_stop('lim_thd_sal : requestead workspace arrays unavailable.')   ;   RETURN
66      END IF
67      ! Set-up pointers to sub-arrays of workspace arrays
68      ze_init =>  wrk_1d_1 (1:jpij)
69      zhiold  =>  wrk_1d_2 (1:jpij)
70      zsiold  =>  wrk_1d_3 (1:jpij)
71
72      !------------------------------------------------------------------------------|
73      ! 1) Constant salinity, constant in time                                       |
74      !------------------------------------------------------------------------------|
75!!gm comment: if num_sal = 1 s_i_b and sm_i_b can be set to bulk_sal one for all in the initialisation phase !!
76      IF( num_sal == 1 ) THEN
77         !
78         DO jk = 1, nlay_i
79            DO ji = kideb, kiut
80               s_i_b(ji,jk) =  bulk_sal
81            END DO ! ji
82         END DO ! jk
83         !
84         DO ji = kideb, kiut
85            sm_i_b(ji)      =  bulk_sal 
86         END DO ! ji
87         !
88      ENDIF
89
90      !------------------------------------------------------------------------------|
91      !  Module 2 : Constant salinity varying in time                                |
92      !------------------------------------------------------------------------------|
93
94      IF(  num_sal == 2  .OR.  num_sal == 4  ) THEN
95
96         !---------------------------------
97         ! Thickness at previous time step
98         !---------------------------------
99         DO ji = kideb, kiut
100            zhiold(ji) = ht_i_b(ji) - dh_i_bott(ji) - dh_snowice(ji) - dh_i_surf(ji)
101         END DO
102
103         !---------------------
104         ! Global heat content
105         !---------------------
106         ze_init(:)  =  0._wp
107         DO jk = 1, nlay_i
108            DO ji = kideb, kiut
109               ze_init(ji) = ze_init(ji) + q_i_b(ji,jk) * ht_i_b(ji) / nlay_i
110            END DO
111         END DO
112
113         DO ji = kideb, kiut
114            !
115            ! Switches
116            !----------
117            iflush       =         MAX( 0._wp , SIGN( 1.0 , t_su_b(ji) - rtt )        )    ! =1 if summer
118            igravdr      =         MAX( 0._wp , SIGN( 1.0 , t_bo_b(ji) - t_su_b(ji) ) )    ! =1 if t_su < t_bo
119            iaccrbo      =         MAX( 0._wp , SIGN( 1.0 , dh_i_bott(ji) )           )    ! =1 if bottom accretion
120            i_ice_switch = 1._wp - MAX ( 0._wp , SIGN( 1._wp , - ht_i_b(ji) + 1.e-2 ) )
121            isnowic      = 1._wp - MAX ( 0._wp , SIGN( 1._wp , - dh_snowice(ji) ) ) * i_ice_switch   ! =1 if snow ice formation
122
123            !---------------------
124            ! Salinity tendencies
125            !---------------------
126            !                                   ! drainage by gravity drainage
127            dsm_i_gd_1d(ji) = - igravdr * MAX( sm_i_b(ji) - sal_G , 0._wp ) / time_G * rdt_ice 
128            !                                   ! drainage by flushing 
129            dsm_i_fl_1d(ji) = - iflush * MAX( sm_i_b(ji) - sal_F , 0._wp ) / time_F * rdt_ice
130
131            !-----------------
132            ! Update salinity   
133            !-----------------
134            ! only drainage terms ( gravity drainage and flushing )
135            ! snow ice / bottom sources are added in lim_thd_ent to conserve energy
136            zsiold(ji) = sm_i_b(ji)
137            sm_i_b(ji) = sm_i_b(ji) + dsm_i_fl_1d(ji) + dsm_i_gd_1d(ji)
138
139            ! if no ice, salinity = 0.1
140            i_ice_switch = 1._wp - MAX ( 0._wp, SIGN( 1._wp , - ht_i_b(ji) ) )
141            sm_i_b(ji)   = i_ice_switch * sm_i_b(ji) + s_i_min * ( 1._wp - i_ice_switch )
142         END DO ! ji
143
144         ! Salinity profile
145         CALL lim_var_salprof1d( kideb, kiut )
146
147         !----------------------------
148         ! Heat flux - brine drainage
149         !----------------------------
150
151         DO ji = kideb, kiut
152!!gm useless
153            ! iflush  : 1 if summer
154            iflush  =  MAX( 0._wp , SIGN ( 1._wp , t_su_b(ji) - rtt ) ) 
155            ! igravdr : 1 if t_su lt t_bo
156            igravdr =  MAX( 0._wp , SIGN ( 1._wp , t_bo_b(ji) - t_su_b(ji) ) ) 
157            ! iaccrbo : 1 if bottom accretion
158            iaccrbo =  MAX( 0._wp , SIGN ( 1._wp , dh_i_bott(ji) ) )
159!!gm end useless
160            !
161            fhbri_1d(ji) = 0._wp
162         END DO ! ji
163
164         !----------------------------
165         ! Salt flux - brine drainage
166         !----------------------------
167         DO ji = kideb, kiut
168            i_ice_switch = 1._wp - MAX ( 0._wp, SIGN( 1._wp , - ht_i_b(ji) ) )
169            fsbri_1d(ji) = fsbri_1d(ji) - i_ice_switch * rhoic * a_i_b(ji) * ht_i_b(ji)         &
170               &         * ( MAX(dsm_i_gd_1d(ji) + dsm_i_fl_1d(ji), sm_i_b(ji) - zsiold(ji) ) ) / rdt_ice
171            IF( num_sal == 4 ) fsbri_1d(ji) = 0._wp
172         END DO ! ji
173
174         ! Only necessary for conservation check since salinity is modified
175         !--------------------
176         ! Temperature update
177         !--------------------
178         DO jk = 1, nlay_i
179            DO ji = kideb, kiut
180               ztmelts    =  -tmut*s_i_b(ji,jk) + rtt
181               !Conversion q(S,T) -> T (second order equation)
182               zaaa         =  cpic
183               zbbb         =  ( rcp - cpic ) * ( ztmelts - rtt ) + q_i_b(ji,jk) / rhoic - lfus
184               zccc         =  lfus * ( ztmelts - rtt )
185               zdiscrim     =  SQRT(  MAX( zbbb*zbbb - 4.0*zaaa*zccc, 0._wp )  )
186               t_i_b(ji,jk) =  rtt - ( zbbb + zdiscrim ) / ( 2.0 *zaaa )
187            END DO
188         END DO
189         !
190      ENDIF ! num_sal .EQ. 2
191
192      !------------------------------------------------------------------------------|
193      !  Module 3 : Profile of salinity, constant in time                            |
194      !------------------------------------------------------------------------------|
195
196      IF( num_sal == 3 )   CALL lim_var_salprof1d( kideb, kiut )
197
198      !------------------------------------------------------------------------------|
199      !  Module 4 : Constant salinity varying in time                                |
200      !------------------------------------------------------------------------------|
201
202      IF( num_sal == 5 ) THEN      ! Cox and Weeks, 1974
203         !
204         DO ji = kideb, kiut
205            zsold = sm_i_b(ji)
206            IF( ht_i_b(ji) < 0.4 ) THEN
207               sm_i_b(ji) = 14.24 - 19.39 * ht_i_b(ji) 
208            ELSE
209               sm_i_b(ji) =  7.88 - 1.59 * ht_i_b(ji)
210               sm_i_b(ji) = MIN( sm_i_b(ji) , zsold ) 
211            ENDIF
212            IF( ht_i_b(ji) > 3.06918239 ) THEN
213               sm_i_b(ji) = 3._wp
214            ENDIF
215            DO jk = 1, nlay_i
216               s_i_b(ji,jk)   = sm_i_b(ji)
217            END DO
218         END DO
219         !
220      ENDIF ! num_sal
221
222      !------------------------------------------------------------------------------|
223      ! 5) Computation of salt flux due to Bottom growth
224      !------------------------------------------------------------------------------|
225
226      IF ( num_sal == 4 ) THEN
227         DO ji = kideb, kiut
228            zji = MOD( npb(ji) - 1 , jpi ) + 1
229            zjj =    ( npb(ji) - 1 ) / jpi + 1
230            fseqv_1d(ji) = fseqv_1d(ji) + ( sss_m(zji,zjj) - bulk_sal    )               &
231               &                        * rhoic * a_i_b(ji) * MAX( dh_i_bott(ji) , 0.0 ) / rdt_ice
232         END DO
233      ELSE
234         DO ji = kideb, kiut
235            zji = MOD( npb(ji) - 1 , jpi ) + 1
236            zjj =    ( npb(ji) - 1 ) / jpi + 1
237            fseqv_1d(ji) = fseqv_1d(ji) + ( sss_m(zji,zjj) - s_i_new(ji) )               &
238               &                        * rhoic * a_i_b(ji) * MAX( dh_i_bott(ji) , 0.0 ) / rdt_ice
239         END DO
240      ENDIF
241      !
242      IF( wrk_not_released(1, 1,2,3) )   CALL ctl_stop( 'lim_thd_sal : failed to release workspace arrays' )
243      !
244   END SUBROUTINE lim_thd_sal
245
246
247   SUBROUTINE lim_thd_sal_init
248      !!-------------------------------------------------------------------
249      !!                  ***  ROUTINE lim_thd_sal_init  ***
250      !!
251      !! ** Purpose :   initialization of ice salinity parameters
252      !!
253      !! ** Method  :   Read the namicesal namelist and check the parameter
254      !!              values called at the first timestep (nit000)
255      !!
256      !! ** input   :   Namelist namicesal
257      !!-------------------------------------------------------------------
258      NAMELIST/namicesal/ num_sal, bulk_sal, sal_G, time_G, sal_F, time_F,   &
259         &                s_i_max, s_i_min, s_i_0, s_i_1
260      !!-------------------------------------------------------------------
261      !
262      REWIND( numnam_ice )                   ! Read Namelist namicesal
263      READ  ( numnam_ice  , namicesal )
264      !
265      IF(lwp) THEN                           ! control print
266         WRITE(numout,*)
267         WRITE(numout,*) 'lim_thd_sal_init : Ice parameters for salinity '
268         WRITE(numout,*) '~~~~~~~~~~~~~~~~'
269         WRITE(numout,*) ' switch for salinity num_sal        : ', num_sal
270         WRITE(numout,*) ' bulk salinity value if num_sal = 1 : ', bulk_sal
271         WRITE(numout,*) ' restoring salinity for GD          : ', sal_G
272         WRITE(numout,*) ' restoring time for GD              : ', time_G
273         WRITE(numout,*) ' restoring salinity for flushing    : ', sal_F
274         WRITE(numout,*) ' restoring time for flushing        : ', time_F
275         WRITE(numout,*) ' Maximum tolerated ice salinity     : ', s_i_max
276         WRITE(numout,*) ' Minimum tolerated ice salinity     : ', s_i_min
277         WRITE(numout,*) ' 1st salinity for salinity profile  : ', s_i_0
278         WRITE(numout,*) ' 2nd salinity for salinity profile  : ', s_i_1
279      ENDIF
280      !
281   END SUBROUTINE lim_thd_sal_init
282
283#else
284   !!----------------------------------------------------------------------
285   !!   Default option         Dummy Module          No LIM-3 sea-ice model
286   !!----------------------------------------------------------------------
287#endif
288   !!======================================================================
289END MODULE limthd_sal
Note: See TracBrowser for help on using the repository browser.