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

source: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90 @ 5059

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

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