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

source: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90 @ 3523

Last change on this file since 3523 was 3523, checked in by gm, 11 years ago

Branch: dev_r3385_NOCS04_HAMF; #665. LIM3 update: bug correction in limtrp and mass flux added in limthd_lac
  • Property svn:keywords set to Id
File size: 10.9 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   USE wrk_nemo         ! work arrays
26
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC   lim_thd_sal        ! called by limthd module
31   PUBLIC   lim_thd_sal_init   ! called by iceini module
32
33   !!----------------------------------------------------------------------
34   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
35   !! $Id$
36   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
37   !!----------------------------------------------------------------------
38CONTAINS
39
40   SUBROUTINE lim_thd_sal( kideb, kiut )
41      !!-------------------------------------------------------------------
42      !!                ***  ROUTINE lim_thd_sal  ***   
43      !!   
44      !! ** Purpose :   computes new salinities in the ice
45      !!
46      !! ** Method  :  3 possibilities
47      !!               -> num_sal = 1 -> Sice = cst    [ice salinity constant in both time & space]
48      !!               -> num_sal = 2 -> Sice = S(z,t) [Vancoppenolle et al. 2005]
49      !!               -> num_sal = 3 -> Sice = S(z)   [multiyear ice]
50      !!---------------------------------------------------------------------
51      INTEGER, INTENT(in) ::   kideb, kiut   ! thickness category index
52      !
53      INTEGER  ::   ji, jk     ! dummy loop indices
54      REAL(wp) ::   zsold, iflush, iaccrbo, igravdr, isnowic, i_ice_switch,  ztmelts   ! local scalars
55      REAL(wp) ::   zaaa, zbbb, zccc, zdiscrim   ! local scalars
56      REAL(wp), POINTER, DIMENSION(:) ::   ze_init, zhiold, zsiold
57      !!---------------------------------------------------------------------
58
59      CALL wrk_alloc( jpij, ze_init, zhiold, zsiold )
60
61      !------------------------------------------------------------------------------|
62      ! 1) Constant salinity, constant in time                                       |
63      !------------------------------------------------------------------------------|
64!!gm comment: if num_sal = 1 s_i_new, s_i_b and sm_i_b can be set to bulk_sal one for all in the initialisation phase !!
65!!gm           ===>>>   simplification of almost all test on num_sal value
66      IF(  num_sal == 1  ) THEN
67            s_i_b  (kideb:kiut,1:nlay_i) =  bulk_sal
68            sm_i_b (kideb:kiut)          =  bulk_sal 
69            s_i_new(kideb:kiut)          =  bulk_sal
70      ENDIF
71
72      !------------------------------------------------------------------------------|
73      !  Module 2 : Constant salinity varying in time                                |
74      !------------------------------------------------------------------------------|
75
76      IF(  num_sal == 2  ) THEN
77
78         !---------------------------------
79         ! Thickness at previous time step
80         !---------------------------------
81         DO ji = kideb, kiut
82            zhiold(ji) = ht_i_b(ji) - dh_i_bott(ji) - dh_snowice(ji) - dh_i_surf(ji)
83         END DO
84
85         !---------------------
86         ! Global heat content
87         !---------------------
88         ze_init(:)  =  0._wp
89         DO jk = 1, nlay_i
90            DO ji = kideb, kiut
91               ze_init(ji) = ze_init(ji) + q_i_b(ji,jk) * ht_i_b(ji) / nlay_i
92            END DO
93         END DO
94
95         DO ji = kideb, kiut
96            !
97            ! Switches
98            !----------
99            iflush       =         MAX( 0._wp , SIGN( 1.0 , t_su_b(ji) - rtt )        )    ! =1 if summer
100            igravdr      =         MAX( 0._wp , SIGN( 1.0 , t_bo_b(ji) - t_su_b(ji) ) )    ! =1 if t_su < t_bo
101            iaccrbo      =         MAX( 0._wp , SIGN( 1.0 , dh_i_bott(ji) )           )    ! =1 if bottom accretion
102            i_ice_switch = 1._wp - MAX ( 0._wp , SIGN( 1._wp , - ht_i_b(ji) + 1.e-2 ) )
103            isnowic      = 1._wp - MAX ( 0._wp , SIGN( 1._wp , - dh_snowice(ji) ) ) * i_ice_switch   ! =1 if snow ice formation
104
105            !---------------------
106            ! Salinity tendencies
107            !---------------------
108            !                                   ! drainage by gravity drainage
109            dsm_i_gd_1d(ji) = - igravdr * MAX( sm_i_b(ji) - sal_G , 0._wp ) / time_G * rdt_ice 
110            !                                   ! drainage by flushing 
111            dsm_i_fl_1d(ji) = - iflush  * MAX( sm_i_b(ji) - sal_F , 0._wp ) / time_F * rdt_ice
112
113            !-----------------
114            ! Update salinity   
115            !-----------------
116            ! only drainage terms ( gravity drainage and flushing )
117            ! snow ice / bottom sources are added in lim_thd_ent to conserve energy
118            zsiold(ji) = sm_i_b(ji)
119            sm_i_b(ji) = sm_i_b(ji) + dsm_i_fl_1d(ji) + dsm_i_gd_1d(ji)
120
121            ! if no ice, salinity = 0.1
122            i_ice_switch = 1._wp - MAX ( 0._wp, SIGN( 1._wp , - ht_i_b(ji) ) )
123            sm_i_b(ji)   = i_ice_switch * sm_i_b(ji) + s_i_min * ( 1._wp - i_ice_switch )
124         END DO ! ji
125
126         CALL lim_var_salprof1d( kideb, kiut )         ! Salinity profile
127
128
129         !----------------------------
130         ! Heat flux - brine drainage
131         !----------------------------
132
133         DO ji = kideb, kiut
134!!gm useless
135            ! iflush  : 1 if summer
136            iflush  =  MAX( 0._wp , SIGN( 1._wp , t_su_b(ji) - rtt        )  ) 
137            ! igravdr : 1 if t_su lt t_bo
138            igravdr =  MAX( 0._wp , SIGN( 1._wp , t_bo_b(ji) - t_su_b(ji) )  ) 
139            ! iaccrbo : 1 if bottom accretion
140            iaccrbo =  MAX( 0._wp , SIGN( 1._wp , dh_i_bott(ji)           )  )
141!!gm end useless
142            !
143            fhbri_1d(ji) = 0._wp
144         END DO ! ji
145
146         !----------------------------
147         ! Salt flux - brine drainage
148         !----------------------------
149         DO ji = kideb, kiut
150            i_ice_switch = 1._wp - MAX(  0._wp, SIGN( 1._wp , - ht_i_b(ji) )  )
151            fsbri_1d(ji) = fsbri_1d(ji) - i_ice_switch * rhoic * a_i_b(ji) * ht_i_b(ji)         &
152               &         * ( MAX( dsm_i_gd_1d(ji) + dsm_i_fl_1d(ji) , sm_i_b(ji) - zsiold(ji) ) ) * r1_rdtice
153         END DO ! ji
154
155         ! Only necessary for conservation check since salinity is modified
156         !--------------------
157         ! Temperature update
158         !--------------------
159         DO jk = 1, nlay_i
160            DO ji = kideb, kiut
161               ztmelts    =  -tmut*s_i_b(ji,jk) + rtt
162               !Conversion q(S,T) -> T (second order equation)
163               zaaa         =  cpic
164               zbbb         =  ( rcp - cpic ) * ( ztmelts - rtt ) + q_i_b(ji,jk) / rhoic - lfus
165               zccc         =  lfus * ( ztmelts - rtt )
166               zdiscrim     =  SQRT(  MAX( zbbb*zbbb - 4.0*zaaa*zccc, 0._wp )  )
167               t_i_b(ji,jk) =  rtt - ( zbbb + zdiscrim ) / ( 2.0 *zaaa )
168            END DO
169         END DO
170         !
171      ENDIF 
172
173      !------------------------------------------------------------------------------|
174      !  Module 3 : Profile of salinity, constant in time                            |
175      !------------------------------------------------------------------------------|
176
177      IF(  num_sal == 3  )   CALL lim_var_salprof1d( kideb, kiut )
178
179
180      !------------------------------------------------------------------------------|
181      ! 5) Computation of salt flux due to Bottom growth
182      !------------------------------------------------------------------------------|
183      ! note: s_i_new = bulk_sal in constant salinity case
184      DO ji = kideb, kiut
185         fseqv_1d(ji) = fseqv_1d(ji) - s_i_new(ji) * rhoic * a_i_b(ji) * MAX( dh_i_bott(ji) , 0._wp ) * r1_rdtice
186      END DO
187      !
188      CALL wrk_dealloc( jpij, ze_init, zhiold, zsiold )
189      !
190   END SUBROUTINE lim_thd_sal
191
192
193   SUBROUTINE lim_thd_sal_init
194      !!-------------------------------------------------------------------
195      !!                  ***  ROUTINE lim_thd_sal_init  ***
196      !!
197      !! ** Purpose :   initialization of ice salinity parameters
198      !!
199      !! ** Method  :   Read the namicesal namelist and check the parameter
200      !!              values called at the first timestep (nit000)
201      !!
202      !! ** input   :   Namelist namicesal
203      !!-------------------------------------------------------------------
204      NAMELIST/namicesal/ num_sal, bulk_sal, sal_G, time_G, sal_F, time_F,   &
205         &                s_i_max, s_i_min, s_i_0, s_i_1
206      !!-------------------------------------------------------------------
207      !
208      REWIND( numnam_ice )                   ! Read Namelist namicesal
209      READ  ( numnam_ice  , namicesal )
210      !
211      IF(lwp) THEN                           ! control print
212         WRITE(numout,*)
213         WRITE(numout,*) 'lim_thd_sal_init : Ice parameters for salinity '
214         WRITE(numout,*) '~~~~~~~~~~~~~~~~'
215         WRITE(numout,*) ' switch for salinity num_sal        : ', num_sal
216         WRITE(numout,*) ' bulk salinity value if num_sal = 1 : ', bulk_sal
217         WRITE(numout,*) ' restoring salinity for GD          : ', sal_G
218         WRITE(numout,*) ' restoring time for GD              : ', time_G
219         WRITE(numout,*) ' restoring salinity for flushing    : ', sal_F
220         WRITE(numout,*) ' restoring time for flushing        : ', time_F
221         WRITE(numout,*) ' Maximum tolerated ice salinity     : ', s_i_max
222         WRITE(numout,*) ' Minimum tolerated ice salinity     : ', s_i_min
223         WRITE(numout,*) ' 1st salinity for salinity profile  : ', s_i_0
224         WRITE(numout,*) ' 2nd salinity for salinity profile  : ', s_i_1
225      ENDIF
226      !
227   END SUBROUTINE lim_thd_sal_init
228
229#else
230   !!----------------------------------------------------------------------
231   !!   Default option         Dummy Module          No LIM-3 sea-ice model
232   !!----------------------------------------------------------------------
233#endif
234   !!======================================================================
235END MODULE limthd_sal
Note: See TracBrowser for help on using the repository browser.