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

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90 @ 3558

Last change on this file since 3558 was 3558, checked in by rblod, 11 years ago

Fix issues when using key_nosignedzeo, see ticket #996

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