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

source: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90 @ 5076

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

LIM3: important bug fix on heat conservation
File size: 10.7 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   !!            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_update1   : 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 wrk_nemo        ! work arrays
26   USE timing          ! Timing
27   USE limcons         ! conservation tests
28   USE lib_mpp         ! MPP library
29   USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
30   USE in_out_manager  ! I/O manager
31
32   IMPLICIT NONE
33   PRIVATE
34
35   PUBLIC   lim_update1
36
37   !! * Substitutions
38#  include "vectopt_loop_substitute.h90"
39   !!----------------------------------------------------------------------
40   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
41   !! $Id: limupdate.F90 3294 2012-01-28 16:44:18Z rblod $
42   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
43   !!----------------------------------------------------------------------
44CONTAINS
45
46   SUBROUTINE lim_update1( kt )
47      !!-------------------------------------------------------------------
48      !!               ***  ROUTINE lim_update1  ***
49      !!               
50      !! ** Purpose :  Computes update of sea-ice global variables at
51      !!               the end of the dynamics.
52      !!               
53      !!---------------------------------------------------------------------
54      INTEGER, INTENT(in) ::   kt    ! number of iteration
55      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
56      REAL(wp) ::   zsal
57      REAL(wp) ::   zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
58      !!-------------------------------------------------------------------
59      IF( nn_timing == 1 )  CALL timing_start('limupdate1')
60
61      IF( ln_limdyn ) THEN
62
63      IF( kt == nit000 .AND. lwp ) THEN
64         WRITE(numout,*) ' lim_update1 ' 
65         WRITE(numout,*) ' ~~~~~~~~~~~ '
66      ENDIF
67
68      ! conservation test
69      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limupdate1', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
70
71      !-----------------
72      ! zap small values
73      !-----------------
74      CALL lim_var_zapsmall
75
76      CALL lim_var_glo2eqv
77     
78      !----------------------------------------------------
79      ! Rebin categories with thickness out of bounds
80      !----------------------------------------------------
81      IF ( jpl > 1 )   CALL lim_itd_th_reb(1, jpl)
82
83      at_i(:,:) = 0._wp
84      DO jl = 1, jpl
85         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
86      END DO
87
88      !----------------------------------------------------
89      ! ice concentration should not exceed amax
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                  ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl)
97               ENDIF
98            END DO
99         END DO
100      END DO
101
102      at_i(:,:) = 0._wp
103      DO jl = 1, jpl
104         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
105      END DO
106   
107      ! --------------------------------------
108      ! Final thickness distribution rebinning
109      ! --------------------------------------
110      IF ( jpl > 1 ) CALL lim_itd_th_reb(1, jpl)
111
112      !-----------------
113      ! zap small values
114      !-----------------
115      CALL lim_var_zapsmall
116
117      !---------------------
118      ! Ice salinity bounds
119      !---------------------
120      IF (  nn_icesal == 2  ) THEN
121         DO jl = 1, jpl
122            DO jj = 1, jpj 
123               DO ji = 1, jpi
124                  zsal            = smv_i(ji,jj,jl)
125                  smv_i(ji,jj,jl) = sm_i(ji,jj,jl) * v_i(ji,jj,jl)
126                  ! salinity stays in bounds
127                  rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, - v_i(ji,jj,jl) ) )
128                  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) )
129                  ! associated salt flux
130                  sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice
131               END DO
132            END DO
133         END DO
134      ENDIF
135
136      ! conservation test
137      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limupdate1', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
138
139      ! -------------------------------------------------
140      ! Diagnostics
141      ! -------------------------------------------------
142      DO jl  = 1, jpl
143         afx_dyn(:,:) = afx_dyn(:,:) + ( a_i(:,:,jl) - a_i_b(:,:,jl) ) * r1_rdtice
144      END DO
145
146      ! heat content variation (W.m-2)
147      DO jj = 1, jpj
148         DO ji = 1, jpi           
149            diag_heat_dhc(ji,jj) = - ( SUM( e_i(ji,jj,1:nlay_i,:) - e_i_b(ji,jj,1:nlay_i,:) ) +  & 
150               &                       SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) )    &
151               &                     ) * r1_rdtice   
152         END DO
153      END DO
154
155      ! -------------------------------------------------
156      ! control prints
157      ! -------------------------------------------------
158      IF(ln_ctl) THEN   ! Control print
159         CALL prt_ctl_info(' ')
160         CALL prt_ctl_info(' - Cell values : ')
161         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
162         CALL prt_ctl(tab2d_1=e12t       , clinfo1=' lim_update1  : cell area   :')
163         CALL prt_ctl(tab2d_1=at_i       , clinfo1=' lim_update1  : at_i        :')
164         CALL prt_ctl(tab2d_1=vt_i       , clinfo1=' lim_update1  : vt_i        :')
165         CALL prt_ctl(tab2d_1=vt_s       , clinfo1=' lim_update1  : vt_s        :')
166         CALL prt_ctl(tab2d_1=strength   , clinfo1=' lim_update1  : strength    :')
167         CALL prt_ctl(tab2d_1=u_ice      , clinfo1=' lim_update1  : u_ice       :', tab2d_2=v_ice      , clinfo2=' v_ice       :')
168         CALL prt_ctl(tab2d_1=u_ice_b    , clinfo1=' lim_update1  : u_ice_b     :', tab2d_2=v_ice_b    , clinfo2=' v_ice_b     :')
169
170         DO jl = 1, jpl
171            CALL prt_ctl_info(' ')
172            CALL prt_ctl_info(' - Category : ', ivar1=jl)
173            CALL prt_ctl_info('   ~~~~~~~~~~')
174            CALL prt_ctl(tab2d_1=ht_i       (:,:,jl)        , clinfo1= ' lim_update1  : ht_i        : ')
175            CALL prt_ctl(tab2d_1=ht_s       (:,:,jl)        , clinfo1= ' lim_update1  : ht_s        : ')
176            CALL prt_ctl(tab2d_1=t_su       (:,:,jl)        , clinfo1= ' lim_update1  : t_su        : ')
177            CALL prt_ctl(tab2d_1=t_s        (:,:,1,jl)      , clinfo1= ' lim_update1  : t_snow      : ')
178            CALL prt_ctl(tab2d_1=sm_i       (:,:,jl)        , clinfo1= ' lim_update1  : sm_i        : ')
179            CALL prt_ctl(tab2d_1=o_i        (:,:,jl)        , clinfo1= ' lim_update1  : o_i         : ')
180            CALL prt_ctl(tab2d_1=a_i        (:,:,jl)        , clinfo1= ' lim_update1  : a_i         : ')
181            CALL prt_ctl(tab2d_1=a_i_b      (:,:,jl)        , clinfo1= ' lim_update1  : a_i_b       : ')
182            CALL prt_ctl(tab2d_1=v_i        (:,:,jl)        , clinfo1= ' lim_update1  : v_i         : ')
183            CALL prt_ctl(tab2d_1=v_i_b      (:,:,jl)        , clinfo1= ' lim_update1  : v_i_b       : ')
184            CALL prt_ctl(tab2d_1=v_s        (:,:,jl)        , clinfo1= ' lim_update1  : v_s         : ')
185            CALL prt_ctl(tab2d_1=v_s_b      (:,:,jl)        , clinfo1= ' lim_update1  : v_s_b       : ')
186            CALL prt_ctl(tab2d_1=e_i        (:,:,1,jl)      , clinfo1= ' lim_update1  : e_i1        : ')
187            CALL prt_ctl(tab2d_1=e_i_b      (:,:,1,jl)      , clinfo1= ' lim_update1  : e_i1_b      : ')
188            CALL prt_ctl(tab2d_1=e_i        (:,:,2,jl)      , clinfo1= ' lim_update1  : e_i2        : ')
189            CALL prt_ctl(tab2d_1=e_i_b      (:,:,2,jl)      , clinfo1= ' lim_update1  : e_i2_b      : ')
190            CALL prt_ctl(tab2d_1=e_s        (:,:,1,jl)      , clinfo1= ' lim_update1  : e_snow      : ')
191            CALL prt_ctl(tab2d_1=e_s_b      (:,:,1,jl)      , clinfo1= ' lim_update1  : e_snow_b    : ')
192            CALL prt_ctl(tab2d_1=smv_i      (:,:,jl)        , clinfo1= ' lim_update1  : smv_i       : ')
193            CALL prt_ctl(tab2d_1=smv_i_b    (:,:,jl)        , clinfo1= ' lim_update1  : smv_i_b     : ')
194            CALL prt_ctl(tab2d_1=oa_i       (:,:,jl)        , clinfo1= ' lim_update1  : oa_i        : ')
195            CALL prt_ctl(tab2d_1=oa_i_b     (:,:,jl)        , clinfo1= ' lim_update1  : oa_i_b      : ')
196
197            DO jk = 1, nlay_i
198               CALL prt_ctl_info(' - Layer : ', ivar1=jk)
199               CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_update1  : t_i       : ')
200            END DO
201         END DO
202
203         CALL prt_ctl_info(' ')
204         CALL prt_ctl_info(' - Heat / FW fluxes : ')
205         CALL prt_ctl_info('   ~~~~~~~~~~~~~~~~~~ ')
206         CALL prt_ctl(tab2d_1=sst_m  , clinfo1= ' lim_update1 : sst   : ', tab2d_2=sss_m     , clinfo2= ' sss       : ')
207
208         CALL prt_ctl_info(' ')
209         CALL prt_ctl_info(' - Stresses : ')
210         CALL prt_ctl_info('   ~~~~~~~~~~ ')
211         CALL prt_ctl(tab2d_1=utau       , clinfo1= ' lim_update1 : utau      : ', tab2d_2=vtau       , clinfo2= ' vtau      : ')
212         CALL prt_ctl(tab2d_1=utau_ice   , clinfo1= ' lim_update1 : utau_ice  : ', tab2d_2=vtau_ice   , clinfo2= ' vtau_ice  : ')
213         CALL prt_ctl(tab2d_1=u_oce      , clinfo1= ' lim_update1 : u_oce     : ', tab2d_2=v_oce      , clinfo2= ' v_oce     : ')
214      ENDIF
215   
216      ENDIF ! ln_limdyn
217
218      IF( nn_timing == 1 )  CALL timing_stop('limupdate1')
219   END SUBROUTINE lim_update1
220#else
221   !!----------------------------------------------------------------------
222   !!   Default option         Empty Module               No sea-ice model
223   !!----------------------------------------------------------------------
224CONTAINS
225   SUBROUTINE lim_update1     ! Empty routine
226   END SUBROUTINE lim_update1
227
228#endif
229
230END MODULE limupdate1
Note: See TracBrowser for help on using the repository browser.