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

source: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90 @ 4655

Last change on this file since 4655 was 4649, checked in by clem, 10 years ago

finalizing LIM3 heat budget conservation + multiple minor bugs corrections
File size: 14.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   !!----------------------------------------------------------------------
8#if defined key_lim3
9   !!----------------------------------------------------------------------
10   !!   'key_lim3'                                      LIM3 sea-ice model
11   !!----------------------------------------------------------------------
12   !!    lim_update2   : computes update of sea-ice global variables from trend terms
13   !!----------------------------------------------------------------------
14   USE limrhg          ! ice rheology
15
16   USE dom_oce
17   USE oce             ! dynamics and tracers variables
18   USE in_out_manager
19   USE sbc_oce         ! Surface boundary condition: ocean fields
20   USE sbc_ice         ! Surface boundary condition: ice fields
21   USE dom_ice
22   USE phycst          ! physical constants
23   USE ice
24   USE limdyn
25   USE limtrp
26   USE limthd
27   USE limsbc
28   USE limdiahsb
29   USE limwri
30   USE limrst
31   USE thd_ice         ! LIM thermodynamic sea-ice variables
32   USE par_ice
33   USE limitd_th
34   USE limitd_me
35   USE limvar
36   USE prtctl           ! Print control
37   USE lbclnk           ! lateral boundary condition - MPP exchanges
38   USE wrk_nemo         ! work arrays
39   USE lib_fortran     ! glob_sum
40   USE timing          ! Timing
41   USE limcons        ! conservation tests
42
43   IMPLICIT NONE
44   PRIVATE
45
46   PUBLIC   lim_update2   ! routine called by ice_step
47
48   REAL(wp)  ::   epsi10 = 1.e-10_wp   !    -       -
49   REAL(wp)  ::   epsi20 = 1.e-20_wp   
50     
51   !! * Substitutions
52#  include "vectopt_loop_substitute.h90"
53   !!----------------------------------------------------------------------
54   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
55   !! $Id: limupdate.F90 3294 2012-01-28 16:44:18Z rblod $
56   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
57   !!----------------------------------------------------------------------
58CONTAINS
59
60   SUBROUTINE lim_update2
61      !!-------------------------------------------------------------------
62      !!               ***  ROUTINE lim_update2  ***
63      !!               
64      !! ** Purpose :  Computes update of sea-ice global variables at
65      !!               the end of the time step.
66      !!               Address pathological cases
67      !!               This place is very important
68      !!               
69      !! ** Method  : 
70      !!    Ice speed from ice dynamics
71      !!    Ice thickness, Snow thickness, Temperatures, Lead fraction
72      !!      from advection and ice thermodynamics
73      !!
74      !! ** Action  : -
75      !!---------------------------------------------------------------------
76      INTEGER ::   ji, jj, jk, jl, jm    ! dummy loop indices
77      INTEGER ::   jbnd1, jbnd2
78      INTEGER ::   i_ice_switch
79      REAL(wp) ::   zindb, zindsn, zindic
80      REAL(wp) ::   zh, zdvres, zsal
81
82      REAL(wp) ::   zEs          ! specific enthalpy of snow (J/kg)
83      REAL(wp) ::   zEi          ! specific enthalpy of ice (J/kg)
84      REAL(wp) ::   zEw          ! specific enthalpy of exchanged water (J/kg)
85      REAL(wp) ::   zdE          ! specific enthalpy difference (J/kg)
86      REAL(wp) ::   zfmdt        ! exchange mass flux x time step (J/m2), >0 towards the ocean
87      !
88      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
89      !!-------------------------------------------------------------------
90      IF( nn_timing == 1 )  CALL timing_start('limupdate2')
91
92      ! conservation test
93      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limupdate2', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
94
95      !-----------------
96      ! zap small values
97      !-----------------
98      CALL lim_itd_me_zapsmall
99
100      CALL lim_var_glo2eqv
101
102      !----------------------------------------------------
103      ! Rebin categories with thickness out of bounds
104      !----------------------------------------------------
105      DO jm = 1, jpm
106         jbnd1 = ice_cat_bounds(jm,1)
107         jbnd2 = ice_cat_bounds(jm,2)
108         IF (ice_ncat_types(jm) .GT. 1 )   CALL lim_itd_th_reb(jbnd1, jbnd2, jm)
109      END DO
110
111      !----------------------------------------------------------------------
112      ! Constrain the thickness of the smallest category above hiclim
113      !----------------------------------------------------------------------
114      DO jm = 1, jpm
115         DO jj = 1, jpj 
116            DO ji = 1, jpi
117               jl = ice_cat_bounds(jm,1)
118               IF( v_i(ji,jj,jl) > 0._wp .AND. ht_i(ji,jj,jl) < hiclim ) THEN
119                  zh             = hiclim / ht_i(ji,jj,jl)
120                  ht_s(ji,jj,jl) = ht_s(ji,jj,jl) * zh
121                  ht_i(ji,jj,jl) = ht_i(ji,jj,jl) * zh
122                  a_i (ji,jj,jl) = a_i(ji,jj,jl)  / zh
123               ENDIF
124            END DO !ji
125         END DO !jj
126      END DO !jm
127     
128      !-----------------------------------------------------
129      ! ice concentration should not exceed amax
130      !-----------------------------------------------------
131      at_i(:,:) = 0._wp
132      DO jl = 1, jpl
133         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
134      END DO
135
136      DO jl  = 1, jpl
137         DO jj = 1, jpj
138            DO ji = 1, jpi
139               IF( at_i(ji,jj) > amax .AND. a_i(ji,jj,jl) > 0._wp ) THEN
140                  a_i(ji,jj,jl)  = a_i(ji,jj,jl) * ( 1._wp - ( 1._wp - amax / at_i(ji,jj) ) )
141                  ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl)
142               ENDIF
143            END DO
144         END DO
145      END DO
146
147      at_i(:,:) = 0.0
148      DO jl = 1, jpl
149         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
150      END DO
151
152      ! --------------------------------------
153      ! Final thickness distribution rebinning
154      ! --------------------------------------
155      DO jm = 1, jpm
156         jbnd1 = ice_cat_bounds(jm,1)
157         jbnd2 = ice_cat_bounds(jm,2)
158         IF (ice_ncat_types(jm) .GT. 1 ) CALL lim_itd_th_reb(jbnd1, jbnd2, jm)
159         IF (ice_ncat_types(jm) .EQ. 1 ) THEN
160         ENDIF
161      END DO
162
163      !-----------------
164      ! zap small values
165      !-----------------
166      CALL lim_itd_me_zapsmall
167
168      !---------------------
169      ! 2.11) Ice salinity
170      !---------------------
171      IF (  num_sal == 2  ) THEN
172         DO jl = 1, jpl
173            DO jj = 1, jpj 
174               DO ji = 1, jpi
175                  zsal            = smv_i(ji,jj,jl)
176                  smv_i(ji,jj,jl) = sm_i(ji,jj,jl) * v_i(ji,jj,jl)
177                  ! salinity stays in bounds
178                  i_ice_switch    = 1._wp - MAX( 0._wp, SIGN( 1._wp, - v_i(ji,jj,jl) ) )
179                  smv_i(ji,jj,jl) = i_ice_switch * 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 - i_ice_switch ) * v_i(ji,jj,jl)
180                  ! associated salt flux
181                  sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice
182               END DO ! ji
183            END DO ! jj
184         END DO !jl
185      ENDIF
186
187      !------------------------------------------------------------------------------
188      ! 2) Corrections to avoid wrong values                                        |
189      !------------------------------------------------------------------------------
190      ! Ice drift
191      !------------
192      DO jj = 2, jpjm1
193         DO ji = 2, jpim1
194            IF ( at_i(ji,jj) == 0._wp ) THEN ! what to do if there is no ice
195               IF ( at_i(ji+1,jj) == 0._wp ) u_ice(ji,jj)   = 0._wp ! right side
196               IF ( at_i(ji-1,jj) == 0._wp ) u_ice(ji-1,jj) = 0._wp ! left side
197               IF ( at_i(ji,jj+1) == 0._wp ) v_ice(ji,jj)   = 0._wp ! upper side
198               IF ( at_i(ji,jj-1) == 0._wp ) v_ice(ji,jj-1) = 0._wp ! bottom side
199            ENDIF
200         END DO
201      END DO
202      !lateral boundary conditions
203      CALL lbc_lnk( u_ice(:,:), 'U', -1. )
204      CALL lbc_lnk( v_ice(:,:), 'V', -1. )
205      !mask velocities
206      u_ice(:,:) = u_ice(:,:) * tmu(:,:)
207      v_ice(:,:) = v_ice(:,:) * tmv(:,:)
208 
209      ! -------------------------------------------------
210      ! Diagnostics
211      ! -------------------------------------------------
212      d_a_i_thd(:,:,:)   = a_i(:,:,:)   - old_a_i(:,:,:) 
213      d_v_s_thd(:,:,:)   = v_s(:,:,:)   - old_v_s(:,:,:)
214      d_v_i_thd(:,:,:)   = v_i(:,:,:)   - old_v_i(:,:,:) 
215      d_e_s_thd(:,:,:,:) = e_s(:,:,:,:) - old_e_s(:,:,:,:) 
216      d_e_i_thd(:,:,1:nlay_i,:) = e_i(:,:,1:nlay_i,:) - old_e_i(:,:,1:nlay_i,:)
217      !?? d_oa_i_thd(:,:,:)  = oa_i (:,:,:) - old_oa_i (:,:,:)
218      d_smv_i_thd(:,:,:) = 0._wp
219      IF( num_sal == 2 )   d_smv_i_thd(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:)
220      ! diag only (clem)
221      dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) * r1_rdtice * rday
222
223      ! heat content variation (W.m-2)
224      DO jj = 1, jpj
225         DO ji = 1, jpi           
226            diag_heat_dhc(ji,jj) = ( SUM( d_e_i_trp(ji,jj,1:nlay_i,:) + d_e_i_thd(ji,jj,1:nlay_i,:) ) +  & 
227               &                     SUM( d_e_s_trp(ji,jj,1:nlay_s,:) + d_e_s_thd(ji,jj,1:nlay_s,:) ) ) * unit_fac * r1_rdtice / area(ji,jj)   
228         END DO
229      END DO
230
231      ! conservation test
232      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limupdate2', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
233
234      IF(ln_ctl) THEN   ! Control print
235         CALL prt_ctl_info(' ')
236         CALL prt_ctl_info(' - Cell values : ')
237         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
238         CALL prt_ctl(tab2d_1=area       , clinfo1=' lim_update2  : cell area   :')
239         CALL prt_ctl(tab2d_1=at_i       , clinfo1=' lim_update2  : at_i        :')
240         CALL prt_ctl(tab2d_1=vt_i       , clinfo1=' lim_update2  : vt_i        :')
241         CALL prt_ctl(tab2d_1=vt_s       , clinfo1=' lim_update2  : vt_s        :')
242         CALL prt_ctl(tab2d_1=strength   , clinfo1=' lim_update2  : strength    :')
243         CALL prt_ctl(tab2d_1=u_ice      , clinfo1=' lim_update2  : u_ice       :', tab2d_2=v_ice      , clinfo2=' v_ice       :')
244         CALL prt_ctl(tab2d_1=old_u_ice  , clinfo1=' lim_update2  : old_u_ice   :', tab2d_2=old_v_ice  , clinfo2=' old_v_ice   :')
245
246         DO jl = 1, jpl
247            CALL prt_ctl_info(' ')
248            CALL prt_ctl_info(' - Category : ', ivar1=jl)
249            CALL prt_ctl_info('   ~~~~~~~~~~')
250            CALL prt_ctl(tab2d_1=ht_i       (:,:,jl)        , clinfo1= ' lim_update2  : ht_i        : ')
251            CALL prt_ctl(tab2d_1=ht_s       (:,:,jl)        , clinfo1= ' lim_update2  : ht_s        : ')
252            CALL prt_ctl(tab2d_1=t_su       (:,:,jl)        , clinfo1= ' lim_update2  : t_su        : ')
253            CALL prt_ctl(tab2d_1=t_s        (:,:,1,jl)      , clinfo1= ' lim_update2  : t_snow      : ')
254            CALL prt_ctl(tab2d_1=sm_i       (:,:,jl)        , clinfo1= ' lim_update2  : sm_i        : ')
255            CALL prt_ctl(tab2d_1=o_i        (:,:,jl)        , clinfo1= ' lim_update2  : o_i         : ')
256            CALL prt_ctl(tab2d_1=a_i        (:,:,jl)        , clinfo1= ' lim_update2  : a_i         : ')
257            CALL prt_ctl(tab2d_1=old_a_i    (:,:,jl)        , clinfo1= ' lim_update2  : old_a_i     : ')
258            CALL prt_ctl(tab2d_1=d_a_i_thd  (:,:,jl)        , clinfo1= ' lim_update2  : d_a_i_thd   : ')
259            CALL prt_ctl(tab2d_1=v_i        (:,:,jl)        , clinfo1= ' lim_update2  : v_i         : ')
260            CALL prt_ctl(tab2d_1=old_v_i    (:,:,jl)        , clinfo1= ' lim_update2  : old_v_i     : ')
261            CALL prt_ctl(tab2d_1=d_v_i_thd  (:,:,jl)        , clinfo1= ' lim_update2  : d_v_i_thd   : ')
262            CALL prt_ctl(tab2d_1=v_s        (:,:,jl)        , clinfo1= ' lim_update2  : v_s         : ')
263            CALL prt_ctl(tab2d_1=old_v_s    (:,:,jl)        , clinfo1= ' lim_update2  : old_v_s     : ')
264            CALL prt_ctl(tab2d_1=d_v_s_thd  (:,:,jl)        , clinfo1= ' lim_update2  : d_v_s_thd   : ')
265            CALL prt_ctl(tab2d_1=e_i        (:,:,1,jl)/1.0e9, clinfo1= ' lim_update2  : e_i1        : ')
266            CALL prt_ctl(tab2d_1=old_e_i    (:,:,1,jl)/1.0e9, clinfo1= ' lim_update2  : old_e_i1    : ')
267            CALL prt_ctl(tab2d_1=d_e_i_thd  (:,:,1,jl)/1.0e9, clinfo1= ' lim_update2  : de_i1_thd   : ')
268            CALL prt_ctl(tab2d_1=e_i        (:,:,2,jl)/1.0e9, clinfo1= ' lim_update2  : e_i2        : ')
269            CALL prt_ctl(tab2d_1=old_e_i    (:,:,2,jl)/1.0e9, clinfo1= ' lim_update2  : old_e_i2    : ')
270            CALL prt_ctl(tab2d_1=d_e_i_thd  (:,:,2,jl)/1.0e9, clinfo1= ' lim_update2  : de_i2_thd   : ')
271            CALL prt_ctl(tab2d_1=e_s        (:,:,1,jl)      , clinfo1= ' lim_update2  : e_snow      : ')
272            CALL prt_ctl(tab2d_1=old_e_s    (:,:,1,jl)      , clinfo1= ' lim_update2  : old_e_snow  : ')
273            CALL prt_ctl(tab2d_1=d_e_s_thd  (:,:,1,jl)/1.0e9, clinfo1= ' lim_update2  : d_e_s_thd   : ')
274            CALL prt_ctl(tab2d_1=smv_i      (:,:,jl)        , clinfo1= ' lim_update2  : smv_i       : ')
275            CALL prt_ctl(tab2d_1=old_smv_i  (:,:,jl)        , clinfo1= ' lim_update2  : old_smv_i   : ')
276            CALL prt_ctl(tab2d_1=d_smv_i_thd(:,:,jl)        , clinfo1= ' lim_update2  : d_smv_i_thd : ')
277            CALL prt_ctl(tab2d_1=oa_i       (:,:,jl)        , clinfo1= ' lim_update2  : oa_i        : ')
278            CALL prt_ctl(tab2d_1=old_oa_i   (:,:,jl)        , clinfo1= ' lim_update2  : old_oa_i    : ')
279            CALL prt_ctl(tab2d_1=d_oa_i_thd (:,:,jl)        , clinfo1= ' lim_update2  : d_oa_i_thd  : ')
280
281            DO jk = 1, nlay_i
282               CALL prt_ctl_info(' - Layer : ', ivar1=jk)
283               CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_update2  : t_i       : ')
284            END DO
285         END DO
286
287         CALL prt_ctl_info(' ')
288         CALL prt_ctl_info(' - Heat / FW fluxes : ')
289         CALL prt_ctl_info('   ~~~~~~~~~~~~~~~~~~ ')
290         CALL prt_ctl(tab2d_1=sst_m  , clinfo1= ' lim_update2 : sst   : ', tab2d_2=sss_m     , clinfo2= ' sss       : ')
291
292         CALL prt_ctl_info(' ')
293         CALL prt_ctl_info(' - Stresses : ')
294         CALL prt_ctl_info('   ~~~~~~~~~~ ')
295         CALL prt_ctl(tab2d_1=utau       , clinfo1= ' lim_update2 : utau      : ', tab2d_2=vtau       , clinfo2= ' vtau      : ')
296         CALL prt_ctl(tab2d_1=utau_ice   , clinfo1= ' lim_update2 : utau_ice  : ', tab2d_2=vtau_ice   , clinfo2= ' vtau_ice  : ')
297         CALL prt_ctl(tab2d_1=u_oce      , clinfo1= ' lim_update2 : u_oce     : ', tab2d_2=v_oce      , clinfo2= ' v_oce     : ')
298      ENDIF
299   
300      IF( nn_timing == 1 )  CALL timing_stop('limupdate2')
301
302   END SUBROUTINE lim_update2
303#else
304   !!----------------------------------------------------------------------
305   !!   Default option         Empty Module               No sea-ice model
306   !!----------------------------------------------------------------------
307CONTAINS
308   SUBROUTINE lim_update2     ! Empty routine
309   END SUBROUTINE lim_update2
310
311#endif
312
313END MODULE limupdate2
Note: See TracBrowser for help on using the repository browser.