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.
limupdate2.F90 in trunk/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90 @ 5134

Last change on this file since 5134 was 5134, checked in by clem, 9 years ago

LIM3: changes to constrain ice thickness after advection. Ice volume and concentration are advected while ice thickness is deduced. This could result in overly high thicknesses associated with very low concentrations (in the 5th category)
File size: 12.1 KB
Line 
1MODULE limupdate2
2   !!======================================================================
3   !!                     ***  MODULE  limupdate2  ***
4   !!   LIM-3 : Update of sea-ice global variables at the end of the time step
5   !!======================================================================
6   !! History :  3.0  !  2006-04  (M. Vancoppenolle) Original code
7   !!            3.5  !  2014-06  (C. Rousset)       Complete rewriting/cleaning
8   !!----------------------------------------------------------------------
9#if defined key_lim3
10   !!----------------------------------------------------------------------
11   !!   'key_lim3'                                      LIM3 sea-ice model
12   !!----------------------------------------------------------------------
13   !!    lim_update2   : computes update of sea-ice global variables from trend terms
14   !!----------------------------------------------------------------------
15   USE sbc_oce         ! Surface boundary condition: ocean fields
16   USE sbc_ice         ! Surface boundary condition: ice fields
17   USE dom_ice
18   USE dom_oce
19   USE phycst          ! physical constants
20   USE ice
21   USE thd_ice         ! LIM thermodynamic sea-ice variables
22   USE limitd_th
23   USE limvar
24   USE prtctl          ! Print control
25   USE lbclnk          ! lateral boundary condition - MPP exchanges
26   USE wrk_nemo        ! work arrays
27   USE timing          ! Timing
28   USE limcons         ! conservation tests
29   USE limctl
30   USE lib_mpp         ! MPP library
31   USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
32   USE in_out_manager
33
34   IMPLICIT NONE
35   PRIVATE
36
37   PUBLIC   lim_update2   ! routine called by ice_step
38
39   !! * Substitutions
40#  include "vectopt_loop_substitute.h90"
41   !!----------------------------------------------------------------------
42   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
43   !! $Id: limupdate.F90 3294 2012-01-28 16:44:18Z rblod $
44   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
45   !!----------------------------------------------------------------------
46CONTAINS
47
48   SUBROUTINE lim_update2( kt )
49      !!-------------------------------------------------------------------
50      !!               ***  ROUTINE lim_update2  ***
51      !!               
52      !! ** Purpose :  Computes update of sea-ice global variables at
53      !!               the end of the time step.
54      !!
55      !!---------------------------------------------------------------------
56      INTEGER, INTENT(in) ::   kt    ! number of iteration
57      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
58      REAL(wp) ::   zsal
59      REAL(wp) ::   zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
60      !!-------------------------------------------------------------------
61      IF( nn_timing == 1 )  CALL timing_start('limupdate2')
62
63      IF( kt == nit000 .AND. lwp ) THEN
64         WRITE(numout,*) ' lim_update2 '
65         WRITE(numout,*) ' ~~~~~~~~~~~ '
66      ENDIF
67
68      ! conservation test
69      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limupdate2', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
70
71      !----------------------------------------------------------------------
72      ! Constrain the thickness of the smallest category above himin
73      !----------------------------------------------------------------------
74      CALL lim_var_glo2eqv
75      DO jj = 1, jpj 
76         DO ji = 1, jpi
77            IF( v_i(ji,jj,1) > 0._wp .AND. ht_i(ji,jj,1) < rn_himin ) THEN
78               a_i (ji,jj,1) = a_i(ji,jj,1) * ht_i(ji,jj,1) / rn_himin
79            ENDIF
80         END DO
81      END DO
82     
83      !-----------------------------------------------------
84      ! ice concentration should not exceed amax
85      !-----------------------------------------------------
86      at_i(:,:) = 0._wp
87      DO jl = 1, jpl
88         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
89      END DO
90
91      DO jl  = 1, jpl
92         DO jj = 1, jpj
93            DO ji = 1, jpi
94               IF( at_i(ji,jj) > rn_amax .AND. a_i(ji,jj,jl) > 0._wp ) THEN
95                  a_i(ji,jj,jl)  = a_i(ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax / at_i(ji,jj) ) )
96               ENDIF
97            END DO
98         END DO
99      END DO
100
101      !----------------------------------------------------
102      ! Rebin categories with thickness out of bounds
103      !----------------------------------------------------
104      IF ( jpl > 1 ) CALL lim_itd_th_reb( 1, jpl )
105
106      !-----------------
107      ! zap small values
108      !-----------------
109      CALL lim_var_zapsmall
110
111      !---------------------
112      ! Ice salinity
113      !---------------------
114      IF (  nn_icesal == 2  ) THEN
115         DO jl = 1, jpl
116            DO jj = 1, jpj 
117               DO ji = 1, jpi
118                  zsal            = smv_i(ji,jj,jl)
119                  smv_i(ji,jj,jl) = sm_i(ji,jj,jl) * v_i(ji,jj,jl)
120                  ! salinity stays in bounds
121                  rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, - v_i(ji,jj,jl) ) )
122                  smv_i(ji,jj,jl) = rswitch * MAX( MIN( rn_simax * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), rn_simin * v_i(ji,jj,jl) )
123                  ! associated salt flux
124                  sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice
125               END DO
126            END DO
127         END DO
128      ENDIF
129
130      !------------------------------------------------------------------------------
131      ! Corrections to avoid wrong values                                        |
132      !------------------------------------------------------------------------------
133      ! Ice drift
134      !------------
135      DO jj = 2, jpjm1
136         DO ji = 2, jpim1
137            IF ( at_i(ji,jj) == 0._wp ) THEN ! what to do if there is no ice
138               IF ( at_i(ji+1,jj) == 0._wp ) u_ice(ji,jj)   = 0._wp ! right side
139               IF ( at_i(ji-1,jj) == 0._wp ) u_ice(ji-1,jj) = 0._wp ! left side
140               IF ( at_i(ji,jj+1) == 0._wp ) v_ice(ji,jj)   = 0._wp ! upper side
141               IF ( at_i(ji,jj-1) == 0._wp ) v_ice(ji,jj-1) = 0._wp ! bottom side
142            ENDIF
143         END DO
144      END DO
145      !lateral boundary conditions
146      CALL lbc_lnk( u_ice(:,:), 'U', -1. )
147      CALL lbc_lnk( v_ice(:,:), 'V', -1. )
148      !mask velocities
149      u_ice(:,:) = u_ice(:,:) * umask(:,:,1)
150      v_ice(:,:) = v_ice(:,:) * vmask(:,:,1)
151 
152      ! for outputs
153      CALL lim_var_glo2eqv            ! equivalent variables (outputs)
154      CALL lim_var_agg(2)             ! aggregate ice thickness categories
155
156      ! conservation test
157      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limupdate2', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
158
159      ! -------------------------------------------------
160      ! Diagnostics
161      ! -------------------------------------------------
162      DO jl  = 1, jpl
163         afx_thd(:,:) = afx_thd(:,:) + ( a_i(:,:,jl) - a_i_b(:,:,jl) ) * r1_rdtice
164      END DO
165      afx_tot = afx_thd + afx_dyn
166
167      ! heat content variation (W.m-2)
168      DO jj = 1, jpj
169         DO ji = 1, jpi           
170            diag_heat_dhc(ji,jj) = diag_heat_dhc(ji,jj) -  &
171               &                   ( SUM( e_i(ji,jj,1:nlay_i,:) - e_i_b(ji,jj,1:nlay_i,:) ) +  & 
172               &                     SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) )    &
173               &                   ) * r1_rdtice   
174         END DO
175      END DO
176
177      ! -------------------------------------------------
178      ! control prints
179      ! -------------------------------------------------
180      IF( ln_icectl )   CALL lim_prt( kt, iiceprt, jiceprt, 2, ' - Final state - ' )   ! control print
181
182      IF(ln_ctl) THEN   ! Control print
183         CALL prt_ctl_info(' ')
184         CALL prt_ctl_info(' - Cell values : ')
185         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
186         CALL prt_ctl(tab2d_1=e12t       , clinfo1=' lim_update2  : cell area   :')
187         CALL prt_ctl(tab2d_1=at_i       , clinfo1=' lim_update2  : at_i        :')
188         CALL prt_ctl(tab2d_1=vt_i       , clinfo1=' lim_update2  : vt_i        :')
189         CALL prt_ctl(tab2d_1=vt_s       , clinfo1=' lim_update2  : vt_s        :')
190         CALL prt_ctl(tab2d_1=strength   , clinfo1=' lim_update2  : strength    :')
191         CALL prt_ctl(tab2d_1=u_ice      , clinfo1=' lim_update2  : u_ice       :', tab2d_2=v_ice      , clinfo2=' v_ice       :')
192         CALL prt_ctl(tab2d_1=u_ice_b    , clinfo1=' lim_update2  : u_ice_b     :', tab2d_2=v_ice_b    , clinfo2=' v_ice_b     :')
193
194         DO jl = 1, jpl
195            CALL prt_ctl_info(' ')
196            CALL prt_ctl_info(' - Category : ', ivar1=jl)
197            CALL prt_ctl_info('   ~~~~~~~~~~')
198            CALL prt_ctl(tab2d_1=ht_i       (:,:,jl)        , clinfo1= ' lim_update2  : ht_i        : ')
199            CALL prt_ctl(tab2d_1=ht_s       (:,:,jl)        , clinfo1= ' lim_update2  : ht_s        : ')
200            CALL prt_ctl(tab2d_1=t_su       (:,:,jl)        , clinfo1= ' lim_update2  : t_su        : ')
201            CALL prt_ctl(tab2d_1=t_s        (:,:,1,jl)      , clinfo1= ' lim_update2  : t_snow      : ')
202            CALL prt_ctl(tab2d_1=sm_i       (:,:,jl)        , clinfo1= ' lim_update2  : sm_i        : ')
203            CALL prt_ctl(tab2d_1=o_i        (:,:,jl)        , clinfo1= ' lim_update2  : o_i         : ')
204            CALL prt_ctl(tab2d_1=a_i        (:,:,jl)        , clinfo1= ' lim_update2  : a_i         : ')
205            CALL prt_ctl(tab2d_1=a_i_b      (:,:,jl)        , clinfo1= ' lim_update2  : a_i_b       : ')
206            CALL prt_ctl(tab2d_1=v_i        (:,:,jl)        , clinfo1= ' lim_update2  : v_i         : ')
207            CALL prt_ctl(tab2d_1=v_i_b      (:,:,jl)        , clinfo1= ' lim_update2  : v_i_b       : ')
208            CALL prt_ctl(tab2d_1=v_s        (:,:,jl)        , clinfo1= ' lim_update2  : v_s         : ')
209            CALL prt_ctl(tab2d_1=v_s_b      (:,:,jl)        , clinfo1= ' lim_update2  : v_s_b       : ')
210            CALL prt_ctl(tab2d_1=e_i        (:,:,1,jl)      , clinfo1= ' lim_update2  : e_i1        : ')
211            CALL prt_ctl(tab2d_1=e_i_b      (:,:,1,jl)      , clinfo1= ' lim_update2  : e_i1_b      : ')
212            CALL prt_ctl(tab2d_1=e_i        (:,:,2,jl)      , clinfo1= ' lim_update2  : e_i2        : ')
213            CALL prt_ctl(tab2d_1=e_i_b      (:,:,2,jl)      , clinfo1= ' lim_update2  : e_i2_b      : ')
214            CALL prt_ctl(tab2d_1=e_s        (:,:,1,jl)      , clinfo1= ' lim_update2  : e_snow      : ')
215            CALL prt_ctl(tab2d_1=e_s_b      (:,:,1,jl)      , clinfo1= ' lim_update2  : e_snow_b    : ')
216            CALL prt_ctl(tab2d_1=smv_i      (:,:,jl)        , clinfo1= ' lim_update2  : smv_i       : ')
217            CALL prt_ctl(tab2d_1=smv_i_b    (:,:,jl)        , clinfo1= ' lim_update2  : smv_i_b     : ')
218            CALL prt_ctl(tab2d_1=oa_i       (:,:,jl)        , clinfo1= ' lim_update2  : oa_i        : ')
219            CALL prt_ctl(tab2d_1=oa_i_b     (:,:,jl)        , clinfo1= ' lim_update2  : oa_i_b      : ')
220
221            DO jk = 1, nlay_i
222               CALL prt_ctl_info(' - Layer : ', ivar1=jk)
223               CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_update2  : t_i       : ')
224            END DO
225         END DO
226
227         CALL prt_ctl_info(' ')
228         CALL prt_ctl_info(' - Heat / FW fluxes : ')
229         CALL prt_ctl_info('   ~~~~~~~~~~~~~~~~~~ ')
230         CALL prt_ctl(tab2d_1=sst_m  , clinfo1= ' lim_update2 : sst   : ', tab2d_2=sss_m     , clinfo2= ' sss       : ')
231
232         CALL prt_ctl_info(' ')
233         CALL prt_ctl_info(' - Stresses : ')
234         CALL prt_ctl_info('   ~~~~~~~~~~ ')
235         CALL prt_ctl(tab2d_1=utau       , clinfo1= ' lim_update2 : utau      : ', tab2d_2=vtau       , clinfo2= ' vtau      : ')
236         CALL prt_ctl(tab2d_1=utau_ice   , clinfo1= ' lim_update2 : utau_ice  : ', tab2d_2=vtau_ice   , clinfo2= ' vtau_ice  : ')
237         CALL prt_ctl(tab2d_1=u_oce      , clinfo1= ' lim_update2 : u_oce     : ', tab2d_2=v_oce      , clinfo2= ' v_oce     : ')
238      ENDIF
239   
240      IF( nn_timing == 1 )  CALL timing_stop('limupdate2')
241
242   END SUBROUTINE lim_update2
243#else
244   !!----------------------------------------------------------------------
245   !!   Default option         Empty Module               No sea-ice model
246   !!----------------------------------------------------------------------
247CONTAINS
248   SUBROUTINE lim_update2     ! Empty routine
249   END SUBROUTINE lim_update2
250
251#endif
252
253END MODULE limupdate2
Note: See TracBrowser for help on using the repository browser.