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.
limupdate1.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/limupdate1.F90 @ 4634

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

major changes in heat budget
File size: 15.9 KB
Line 
1MODULE limupdate1
2   !!======================================================================
3   !!                     ***  MODULE  limupdate1  ***
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_update1   : 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   ! Check budget (Rousset)
41   USE in_out_manager   ! I/O manager
42   USE iom              ! I/O manager
43   USE lib_mpp          ! MPP library
44   USE timing          ! Timing
45
46   IMPLICIT NONE
47   PRIVATE
48
49   PUBLIC   lim_update1   ! routine called by ice_step
50
51      REAL(wp)  ::   epsi10 = 1.e-10_wp   !    -       -
52         
53   !! * Substitutions
54#  include "vectopt_loop_substitute.h90"
55   !!----------------------------------------------------------------------
56   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
57   !! $Id: limupdate.F90 3294 2012-01-28 16:44:18Z rblod $
58   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
59   !!----------------------------------------------------------------------
60CONTAINS
61
62   SUBROUTINE lim_update1
63      !!-------------------------------------------------------------------
64      !!               ***  ROUTINE lim_update1  ***
65      !!               
66      !! ** Purpose :  Computes update of sea-ice global variables at
67      !!               the end of the time step.
68      !!               Address pathological cases
69      !!               This place is very important
70      !!               
71      !! ** Method  : 
72      !!    Ice speed from ice dynamics
73      !!    Ice thickness, Snow thickness, Temperatures, Lead fraction
74      !!      from advection and ice thermodynamics
75      !!
76      !! ** Action  : -
77      !!---------------------------------------------------------------------
78      INTEGER ::   ji, jj, jk, jl, jm    ! dummy loop indices
79      INTEGER ::   jbnd1, jbnd2
80      INTEGER ::   i_ice_switch
81      REAL(wp) ::   zinda, zindb, zindsn, zindic
82      REAL(wp) ::   zh, zdvres, zsal
83      REAL(wp) ::   zat_i_old, ztmelts
84
85      REAL(wp) :: zchk_v_i, zchk_smv, zchk_e_i, zchk_fs, zchk_fw, zchk_ft, zchk_v_i_b, zchk_smv_b, zchk_e_i_b, zchk_fs_b, zchk_fw_b, zchk_ft_b ! Check conservation (C Rousset)
86      REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset)
87      !!-------------------------------------------------------------------
88      IF( nn_timing == 1 )  CALL timing_start('limupdate1')
89
90      !------------------------------------------------------------------------------
91      ! 1. Update of Global variables                                               |
92      !------------------------------------------------------------------------------
93
94      !-----------------
95      !  Trend terms
96      !-----------------
97
98      ! -------------------------------
99      !- check conservation (C Rousset)
100      IF (ln_limdiahsb) THEN
101         zchk_v_i_b = glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) )
102         zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) )
103         zchk_e_i_b = glob_sum( SUM(   e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) )
104         zchk_fw_b  = glob_sum( ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) ) * area(:,:) * tms(:,:) )
105         zchk_fs_b  = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) )
106         zchk_ft_b  = glob_sum( ( hfx_tot(:,:) - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) ) * area(:,:) / unit_fac * tms(:,:) )
107      ENDIF
108      !- check conservation (C Rousset)
109      ! -------------------------------
110
111      ! zap small values
112      !-----------------
113      CALL lim_itd_me_zapsmall
114
115      CALL lim_var_glo2eqv
116     
117      at_i(:,:) = 0._wp
118      DO jl = 1, jpl
119         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
120      END DO
121
122      !----------------------------------------------------
123      ! 2.2) Rebin categories with thickness out of bounds
124      !----------------------------------------------------
125      DO jm = 1, jpm
126         jbnd1 = ice_cat_bounds(jm,1)
127         jbnd2 = ice_cat_bounds(jm,2)
128         IF (ice_ncat_types(jm) .GT. 1 )   CALL lim_itd_th_reb(jbnd1, jbnd2, jm)
129      END DO
130
131      at_i(:,:) = 0._wp
132      DO jl = 1, jpl
133         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
134      END DO
135
136      !--- 2.13 ice concentration should not exceed amax
137      !-----------------------------------------------------
138      DO jl  = 1, jpl
139         DO jj = 1, jpj
140            DO ji = 1, jpi
141               IF( at_i(ji,jj) > amax .AND. a_i(ji,jj,jl) > 0._wp ) THEN
142                  a_i(ji,jj,jl)  = a_i(ji,jj,jl) * ( 1._wp - ( 1._wp - amax / at_i(ji,jj) ) )
143                  ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl)
144               ENDIF
145            END DO
146         END DO
147      END DO
148
149      at_i(:,:) = 0.0
150      DO jl = 1, jpl
151         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
152      END DO
153 
154 
155      ! Final thickness distribution rebinning
156      ! --------------------------------------
157      DO jm = 1, jpm
158         jbnd1 = ice_cat_bounds(jm,1)
159         jbnd2 = ice_cat_bounds(jm,2)
160         IF (ice_ncat_types(jm) .GT. 1 ) CALL lim_itd_th_reb(jbnd1, jbnd2, jm)
161         IF (ice_ncat_types(jm) .EQ. 1 ) THEN
162         ENDIF
163      END DO
164
165
166      ! zap small values
167      !-----------------
168      CALL lim_itd_me_zapsmall
169
170      !---------------------
171      ! 2.11) Ice salinity
172      !---------------------
173      IF (  num_sal == 2  ) THEN
174         DO jl = 1, jpl
175            DO jj = 1, jpj 
176               DO ji = 1, jpi
177                  zsal            = smv_i(ji,jj,jl)
178                  smv_i(ji,jj,jl) = sm_i(ji,jj,jl) * v_i(ji,jj,jl)
179                  ! salinity stays in bounds
180                  i_ice_switch    = 1._wp - MAX( 0._wp, SIGN( 1._wp, - v_i(ji,jj,jl) ) )
181                  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)
182                  ! associated salt flux
183                  sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice
184               END DO ! ji
185            END DO ! jj
186         END DO !jl
187      ENDIF
188
189      ! -------------------
190      at_i(:,:) = a_i(:,:,1)
191      DO jl = 2, jpl
192         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
193      END DO
194 
195
196      ! -------------------------------------------------
197      ! Diagnostics
198      ! -------------------------------------------------
199      d_u_ice_dyn(:,:)     = u_ice(:,:)     - old_u_ice(:,:)
200      d_v_ice_dyn(:,:)     = v_ice(:,:)     - old_v_ice(:,:)
201      d_a_i_trp  (:,:,:)   = a_i  (:,:,:)   - old_a_i  (:,:,:)
202      d_v_s_trp  (:,:,:)   = v_s  (:,:,:)   - old_v_s  (:,:,:) 
203      d_v_i_trp  (:,:,:)   = v_i  (:,:,:)   - old_v_i  (:,:,:)   
204      d_e_s_trp  (:,:,:,:) = e_s  (:,:,:,:) - old_e_s  (:,:,:,:) 
205      d_e_i_trp  (:,:,1:nlay_i,:) = e_i  (:,:,1:nlay_i,:) - old_e_i(:,:,1:nlay_i,:)
206      d_oa_i_trp (:,:,:)   = oa_i (:,:,:)   - old_oa_i (:,:,:)
207      d_smv_i_trp(:,:,:)   = 0._wp
208      IF(  num_sal == 2  )   d_smv_i_trp(:,:,:)  = smv_i(:,:,:) - old_smv_i(:,:,:)
209
210      ! -------------------------------
211      !- check conservation (C Rousset)
212      IF( ln_limdiahsb ) THEN
213         zchk_fs  = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b
214         zchk_fw  = glob_sum( ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fw_b
215         zchk_ft  = glob_sum( ( hfx_tot(:,:) - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) ) * area(:,:) / unit_fac * tms(:,:) ) - zchk_ft_b
216 
217         zchk_v_i = ( glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b ) * r1_rdtice - zchk_fw 
218         zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic )
219         zchk_e_i =   glob_sum( SUM( e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) ) * r1_rdtice - zchk_e_i_b * r1_rdtice + zchk_ft
220
221         zchk_vmin = glob_min(v_i)
222         zchk_amax = glob_max(SUM(a_i,dim=3))
223         zchk_amin = glob_min(a_i)
224
225         IF(lwp) THEN
226            IF ( ABS( zchk_v_i   ) >  1.e-4 ) WRITE(numout,*) 'violation volume [kg/day]     (limupdate1) = ',(zchk_v_i * rday)
227            IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limupdate1) = ',(zchk_smv * rday)
228            IF ( ABS( zchk_e_i   ) >  1.e-2 ) WRITE(numout,*) 'violation enthalpy [1e9 J]    (limupdate1) = ',(zchk_e_i)
229            IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limupdate1) = ',(zchk_vmin * 1.e-3)
230            IF ( zchk_amax >  amax+epsi10   ) WRITE(numout,*) 'violation a_i>amax            (limupdate1) = ',zchk_amax
231            IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limupdate1) = ',zchk_amin
232         ENDIF
233       ENDIF
234      !- check conservation (C Rousset)
235      ! -------------------------------
236
237      IF(ln_ctl) THEN   ! Control print
238         CALL prt_ctl_info(' ')
239         CALL prt_ctl_info(' - Cell values : ')
240         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
241         CALL prt_ctl(tab2d_1=area       , clinfo1=' lim_update1  : cell area   :')
242         CALL prt_ctl(tab2d_1=at_i       , clinfo1=' lim_update1  : at_i        :')
243         CALL prt_ctl(tab2d_1=vt_i       , clinfo1=' lim_update1  : vt_i        :')
244         CALL prt_ctl(tab2d_1=vt_s       , clinfo1=' lim_update1  : vt_s        :')
245         CALL prt_ctl(tab2d_1=strength   , clinfo1=' lim_update1  : strength    :')
246         CALL prt_ctl(tab2d_1=u_ice      , clinfo1=' lim_update1  : u_ice       :', tab2d_2=v_ice      , clinfo2=' v_ice       :')
247         CALL prt_ctl(tab2d_1=d_u_ice_dyn, clinfo1=' lim_update1  : d_u_ice_dyn :', tab2d_2=d_v_ice_dyn, clinfo2=' d_v_ice_dyn :')
248         CALL prt_ctl(tab2d_1=old_u_ice  , clinfo1=' lim_update1  : old_u_ice   :', tab2d_2=old_v_ice  , clinfo2=' old_v_ice   :')
249
250         DO jl = 1, jpl
251            CALL prt_ctl_info(' ')
252            CALL prt_ctl_info(' - Category : ', ivar1=jl)
253            CALL prt_ctl_info('   ~~~~~~~~~~')
254            CALL prt_ctl(tab2d_1=ht_i       (:,:,jl)        , clinfo1= ' lim_update1  : ht_i        : ')
255            CALL prt_ctl(tab2d_1=ht_s       (:,:,jl)        , clinfo1= ' lim_update1  : ht_s        : ')
256            CALL prt_ctl(tab2d_1=t_su       (:,:,jl)        , clinfo1= ' lim_update1  : t_su        : ')
257            CALL prt_ctl(tab2d_1=t_s        (:,:,1,jl)      , clinfo1= ' lim_update1  : t_snow      : ')
258            CALL prt_ctl(tab2d_1=sm_i       (:,:,jl)        , clinfo1= ' lim_update1  : sm_i        : ')
259            CALL prt_ctl(tab2d_1=o_i        (:,:,jl)        , clinfo1= ' lim_update1  : o_i         : ')
260            CALL prt_ctl(tab2d_1=a_i        (:,:,jl)        , clinfo1= ' lim_update1  : a_i         : ')
261            CALL prt_ctl(tab2d_1=old_a_i    (:,:,jl)        , clinfo1= ' lim_update1  : old_a_i     : ')
262            CALL prt_ctl(tab2d_1=d_a_i_trp  (:,:,jl)        , clinfo1= ' lim_update1  : d_a_i_trp   : ')
263            CALL prt_ctl(tab2d_1=v_i        (:,:,jl)        , clinfo1= ' lim_update1  : v_i         : ')
264            CALL prt_ctl(tab2d_1=old_v_i    (:,:,jl)        , clinfo1= ' lim_update1  : old_v_i     : ')
265            CALL prt_ctl(tab2d_1=d_v_i_trp  (:,:,jl)        , clinfo1= ' lim_update1  : d_v_i_trp   : ')
266            CALL prt_ctl(tab2d_1=v_s        (:,:,jl)        , clinfo1= ' lim_update1  : v_s         : ')
267            CALL prt_ctl(tab2d_1=old_v_s    (:,:,jl)        , clinfo1= ' lim_update1  : old_v_s     : ')
268            CALL prt_ctl(tab2d_1=d_v_s_trp  (:,:,jl)        , clinfo1= ' lim_update1  : d_v_s_trp   : ')
269            CALL prt_ctl(tab2d_1=e_i        (:,:,1,jl)/1.0e9, clinfo1= ' lim_update1  : e_i1        : ')
270            CALL prt_ctl(tab2d_1=old_e_i    (:,:,1,jl)/1.0e9, clinfo1= ' lim_update1  : old_e_i1    : ')
271            CALL prt_ctl(tab2d_1=d_e_i_trp  (:,:,1,jl)/1.0e9, clinfo1= ' lim_update1  : de_i1_trp   : ')
272            CALL prt_ctl(tab2d_1=e_i        (:,:,2,jl)/1.0e9, clinfo1= ' lim_update1  : e_i2        : ')
273            CALL prt_ctl(tab2d_1=old_e_i    (:,:,2,jl)/1.0e9, clinfo1= ' lim_update1  : old_e_i2    : ')
274            CALL prt_ctl(tab2d_1=d_e_i_trp  (:,:,2,jl)/1.0e9, clinfo1= ' lim_update1  : de_i2_trp   : ')
275            CALL prt_ctl(tab2d_1=e_s        (:,:,1,jl)      , clinfo1= ' lim_update1  : e_snow      : ')
276            CALL prt_ctl(tab2d_1=old_e_s    (:,:,1,jl)      , clinfo1= ' lim_update1  : old_e_snow  : ')
277            CALL prt_ctl(tab2d_1=d_e_s_trp  (:,:,1,jl)/1.0e9, clinfo1= ' lim_update1  : d_e_s_trp   : ')
278            CALL prt_ctl(tab2d_1=smv_i      (:,:,jl)        , clinfo1= ' lim_update1  : smv_i       : ')
279            CALL prt_ctl(tab2d_1=old_smv_i  (:,:,jl)        , clinfo1= ' lim_update1  : old_smv_i   : ')
280            CALL prt_ctl(tab2d_1=d_smv_i_trp(:,:,jl)        , clinfo1= ' lim_update1  : d_smv_i_trp : ')
281            CALL prt_ctl(tab2d_1=oa_i       (:,:,jl)        , clinfo1= ' lim_update1  : oa_i        : ')
282            CALL prt_ctl(tab2d_1=old_oa_i   (:,:,jl)        , clinfo1= ' lim_update1  : old_oa_i    : ')
283            CALL prt_ctl(tab2d_1=d_oa_i_trp (:,:,jl)        , clinfo1= ' lim_update1  : d_oa_i_trp  : ')
284
285            DO jk = 1, nlay_i
286               CALL prt_ctl_info(' - Layer : ', ivar1=jk)
287               CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_update1  : t_i       : ')
288            END DO
289         END DO
290
291         CALL prt_ctl_info(' ')
292         CALL prt_ctl_info(' - Heat / FW fluxes : ')
293         CALL prt_ctl_info('   ~~~~~~~~~~~~~~~~~~ ')
294         CALL prt_ctl(tab2d_1=sst_m  , clinfo1= ' lim_update1 : sst   : ', tab2d_2=sss_m     , clinfo2= ' sss       : ')
295
296         CALL prt_ctl_info(' ')
297         CALL prt_ctl_info(' - Stresses : ')
298         CALL prt_ctl_info('   ~~~~~~~~~~ ')
299         CALL prt_ctl(tab2d_1=utau       , clinfo1= ' lim_update1 : utau      : ', tab2d_2=vtau       , clinfo2= ' vtau      : ')
300         CALL prt_ctl(tab2d_1=utau_ice   , clinfo1= ' lim_update1 : utau_ice  : ', tab2d_2=vtau_ice   , clinfo2= ' vtau_ice  : ')
301         CALL prt_ctl(tab2d_1=u_oce      , clinfo1= ' lim_update1 : u_oce     : ', tab2d_2=v_oce      , clinfo2= ' v_oce     : ')
302      ENDIF
303   
304      IF( nn_timing == 1 )  CALL timing_stop('limupdate1')
305   END SUBROUTINE lim_update1
306#else
307   !!----------------------------------------------------------------------
308   !!   Default option         Empty Module               No sea-ice model
309   !!----------------------------------------------------------------------
310CONTAINS
311   SUBROUTINE lim_update1     ! Empty routine
312   END SUBROUTINE lim_update1
313
314#endif
315
316END MODULE limupdate1
Note: See TracBrowser for help on using the repository browser.