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 NEMO/releases/release-3.6/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: NEMO/releases/release-3.6/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90 @ 10379

Last change on this file since 10379 was 10379, checked in by clem, 6 years ago

bug fix affecting the lead fraction transfered from the ice model to BGC. This fraction appears not to be limited by 1 (or amax to be more precise) while it should be.

  • Property svn:keywords set to Id
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 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$
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,*)''
65         WRITE(numout,*)' lim_update2 '
66         WRITE(numout,*)' ~~~~~~~~~~~ '
67      ENDIF
68
69      ! conservation test
70      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limupdate2', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
71
72      !----------------------------------------------------------------------
73      ! Constrain the thickness of the smallest category above himin
74      !----------------------------------------------------------------------
75      DO jj = 1, jpj 
76         DO ji = 1, jpi
77            rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,1) - epsi20 ) )   !0 if no ice and 1 if yes
78            ht_i(ji,jj,1) = v_i (ji,jj,1) / MAX( a_i(ji,jj,1) , epsi20 ) * rswitch
79            IF( v_i(ji,jj,1) > 0._wp .AND. ht_i(ji,jj,1) < rn_himin ) THEN
80               a_i (ji,jj,1) = a_i (ji,jj,1) * ht_i(ji,jj,1) / rn_himin
81               oa_i(ji,jj,1) = oa_i(ji,jj,1) * ht_i(ji,jj,1) / rn_himin
82            ENDIF
83         END DO
84      END DO
85     
86      !-----------------------------------------------------
87      ! ice concentration should not exceed amax
88      !-----------------------------------------------------
89      at_i(:,:) = 0._wp
90      DO jl = 1, jpl
91         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
92      END DO
93
94      DO jl  = 1, jpl
95         DO jj = 1, jpj
96            DO ji = 1, jpi
97               IF( at_i(ji,jj) > rn_amax_2d(ji,jj) .AND. a_i(ji,jj,jl) > 0._wp ) THEN
98                  a_i (ji,jj,jl) = a_i (ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax_2d(ji,jj) / at_i(ji,jj) ) )
99                  oa_i(ji,jj,jl) = oa_i(ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax_2d(ji,jj) / at_i(ji,jj) ) )
100               ENDIF
101            END DO
102         END DO
103      END DO
104      at_i(:,:) = SUM( a_i(:,:,:), dim=3 )
105
106      !---------------------
107      ! Ice salinity
108      !---------------------
109      IF (  nn_icesal == 2  ) THEN
110         DO jl = 1, jpl
111            DO jj = 1, jpj 
112               DO ji = 1, jpi
113                  zsal            = smv_i(ji,jj,jl)
114                  ! salinity stays in bounds
115                  rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, - v_i(ji,jj,jl) ) )
116                  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) )
117                  ! associated salt flux
118                  sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice
119               END DO
120            END DO
121         END DO
122      ENDIF
123
124      !----------------------------------------------------
125      ! Rebin categories with thickness out of bounds
126      !----------------------------------------------------
127      IF ( jpl > 1 ) CALL lim_itd_th_reb( 1, jpl )
128
129      !-----------------
130      ! zap small values
131      !-----------------
132      CALL lim_var_zapsmall
133
134      !------------------------------------------------------------------------------
135      ! Corrections to avoid wrong values                                        |
136      !------------------------------------------------------------------------------
137      ! Ice drift
138      !------------
139      DO jj = 2, jpjm1
140         DO ji = 2, jpim1
141            IF ( at_i(ji,jj) == 0._wp ) THEN ! what to do if there is no ice
142               IF ( at_i(ji+1,jj) == 0._wp ) u_ice(ji,jj)   = 0._wp ! right side
143               IF ( at_i(ji-1,jj) == 0._wp ) u_ice(ji-1,jj) = 0._wp ! left side
144               IF ( at_i(ji,jj+1) == 0._wp ) v_ice(ji,jj)   = 0._wp ! upper side
145               IF ( at_i(ji,jj-1) == 0._wp ) v_ice(ji,jj-1) = 0._wp ! bottom side
146            ENDIF
147         END DO
148      END DO
149      !lateral boundary conditions
150      CALL lbc_lnk( u_ice(:,:), 'U', -1. )
151      CALL lbc_lnk( v_ice(:,:), 'V', -1. )
152      !mask velocities
153      u_ice(:,:) = u_ice(:,:) * umask(:,:,1)
154      v_ice(:,:) = v_ice(:,:) * vmask(:,:,1)
155 
156      ! -------------------------------------------------
157      ! Diagnostics
158      ! -------------------------------------------------
159      DO jl  = 1, jpl
160         oa_i(:,:,jl) = oa_i(:,:,jl) + a_i(:,:,jl) * rdt_ice / rday   ! ice natural aging
161         afx_thd(:,:) = afx_thd(:,:) + ( a_i(:,:,jl) - a_i_b(:,:,jl) ) * r1_rdtice
162      END DO
163      afx_tot = afx_thd + afx_dyn
164
165      DO jj = 1, jpj
166         DO ji = 1, jpi           
167            ! heat content variation (W.m-2)
168            diag_heat(ji,jj) = diag_heat(ji,jj) -  &
169               &               ( SUM( e_i(ji,jj,1:nlay_i,:) - e_i_b(ji,jj,1:nlay_i,:) ) +  & 
170               &                 SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) )    &
171               &               ) * r1_rdtice   
172            ! salt, volume
173            diag_smvi(ji,jj) = diag_smvi(ji,jj) + SUM( smv_i(ji,jj,:) - smv_i_b(ji,jj,:) ) * rhoic * r1_rdtice
174            diag_vice(ji,jj) = diag_vice(ji,jj) + SUM( v_i  (ji,jj,:) - v_i_b  (ji,jj,:) ) * rhoic * r1_rdtice
175            diag_vsnw(ji,jj) = diag_vsnw(ji,jj) + SUM( v_s  (ji,jj,:) - v_s_b  (ji,jj,:) ) * rhosn * r1_rdtice
176         END DO
177      END DO
178
179      ! conservation test
180      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limupdate2', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
181
182      ! necessary calls (at least for coupling)
183      CALL lim_var_glo2eqv
184      CALL lim_var_agg(2)
185
186      ! -------------------------------------------------
187      ! control prints
188      ! -------------------------------------------------
189      IF( ln_icectl )   CALL lim_prt( kt, iiceprt, jiceprt, 2, ' - Final state - ' )   ! control print
190
191      IF(ln_ctl) THEN   ! Control print
192         CALL prt_ctl_info(' ')
193         CALL prt_ctl_info(' - Cell values : ')
194         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
195         CALL prt_ctl(tab2d_1=e12t       , clinfo1=' lim_update2  : cell area   :')
196         CALL prt_ctl(tab2d_1=at_i       , clinfo1=' lim_update2  : at_i        :')
197         CALL prt_ctl(tab2d_1=vt_i       , clinfo1=' lim_update2  : vt_i        :')
198         CALL prt_ctl(tab2d_1=vt_s       , clinfo1=' lim_update2  : vt_s        :')
199         CALL prt_ctl(tab2d_1=strength   , clinfo1=' lim_update2  : strength    :')
200         CALL prt_ctl(tab2d_1=u_ice      , clinfo1=' lim_update2  : u_ice       :', tab2d_2=v_ice      , clinfo2=' v_ice       :')
201         CALL prt_ctl(tab2d_1=u_ice_b    , clinfo1=' lim_update2  : u_ice_b     :', tab2d_2=v_ice_b    , clinfo2=' v_ice_b     :')
202
203         DO jl = 1, jpl
204            CALL prt_ctl_info(' ')
205            CALL prt_ctl_info(' - Category : ', ivar1=jl)
206            CALL prt_ctl_info('   ~~~~~~~~~~')
207            CALL prt_ctl(tab2d_1=ht_i       (:,:,jl)        , clinfo1= ' lim_update2  : ht_i        : ')
208            CALL prt_ctl(tab2d_1=ht_s       (:,:,jl)        , clinfo1= ' lim_update2  : ht_s        : ')
209            CALL prt_ctl(tab2d_1=t_su       (:,:,jl)        , clinfo1= ' lim_update2  : t_su        : ')
210            CALL prt_ctl(tab2d_1=t_s        (:,:,1,jl)      , clinfo1= ' lim_update2  : t_snow      : ')
211            CALL prt_ctl(tab2d_1=sm_i       (:,:,jl)        , clinfo1= ' lim_update2  : sm_i        : ')
212            CALL prt_ctl(tab2d_1=o_i        (:,:,jl)        , clinfo1= ' lim_update2  : o_i         : ')
213            CALL prt_ctl(tab2d_1=a_i        (:,:,jl)        , clinfo1= ' lim_update2  : a_i         : ')
214            CALL prt_ctl(tab2d_1=a_i_b      (:,:,jl)        , clinfo1= ' lim_update2  : a_i_b       : ')
215            CALL prt_ctl(tab2d_1=v_i        (:,:,jl)        , clinfo1= ' lim_update2  : v_i         : ')
216            CALL prt_ctl(tab2d_1=v_i_b      (:,:,jl)        , clinfo1= ' lim_update2  : v_i_b       : ')
217            CALL prt_ctl(tab2d_1=v_s        (:,:,jl)        , clinfo1= ' lim_update2  : v_s         : ')
218            CALL prt_ctl(tab2d_1=v_s_b      (:,:,jl)        , clinfo1= ' lim_update2  : v_s_b       : ')
219            CALL prt_ctl(tab2d_1=e_i        (:,:,1,jl)      , clinfo1= ' lim_update2  : e_i1        : ')
220            CALL prt_ctl(tab2d_1=e_i_b      (:,:,1,jl)      , clinfo1= ' lim_update2  : e_i1_b      : ')
221            CALL prt_ctl(tab2d_1=e_i        (:,:,2,jl)      , clinfo1= ' lim_update2  : e_i2        : ')
222            CALL prt_ctl(tab2d_1=e_i_b      (:,:,2,jl)      , clinfo1= ' lim_update2  : e_i2_b      : ')
223            CALL prt_ctl(tab2d_1=e_s        (:,:,1,jl)      , clinfo1= ' lim_update2  : e_snow      : ')
224            CALL prt_ctl(tab2d_1=e_s_b      (:,:,1,jl)      , clinfo1= ' lim_update2  : e_snow_b    : ')
225            CALL prt_ctl(tab2d_1=smv_i      (:,:,jl)        , clinfo1= ' lim_update2  : smv_i       : ')
226            CALL prt_ctl(tab2d_1=smv_i_b    (:,:,jl)        , clinfo1= ' lim_update2  : smv_i_b     : ')
227            CALL prt_ctl(tab2d_1=oa_i       (:,:,jl)        , clinfo1= ' lim_update2  : oa_i        : ')
228            CALL prt_ctl(tab2d_1=oa_i_b     (:,:,jl)        , clinfo1= ' lim_update2  : oa_i_b      : ')
229
230            DO jk = 1, nlay_i
231               CALL prt_ctl_info(' - Layer : ', ivar1=jk)
232               CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_update2  : t_i       : ')
233            END DO
234         END DO
235
236         CALL prt_ctl_info(' ')
237         CALL prt_ctl_info(' - Heat / FW fluxes : ')
238         CALL prt_ctl_info('   ~~~~~~~~~~~~~~~~~~ ')
239         CALL prt_ctl(tab2d_1=sst_m  , clinfo1= ' lim_update2 : sst   : ', tab2d_2=sss_m     , clinfo2= ' sss       : ')
240
241         CALL prt_ctl_info(' ')
242         CALL prt_ctl_info(' - Stresses : ')
243         CALL prt_ctl_info('   ~~~~~~~~~~ ')
244         CALL prt_ctl(tab2d_1=utau       , clinfo1= ' lim_update2 : utau      : ', tab2d_2=vtau       , clinfo2= ' vtau      : ')
245         CALL prt_ctl(tab2d_1=utau_ice   , clinfo1= ' lim_update2 : utau_ice  : ', tab2d_2=vtau_ice   , clinfo2= ' vtau_ice  : ')
246         CALL prt_ctl(tab2d_1=u_oce      , clinfo1= ' lim_update2 : u_oce     : ', tab2d_2=v_oce      , clinfo2= ' v_oce     : ')
247      ENDIF
248   
249      IF( nn_timing == 1 )  CALL timing_stop('limupdate2')
250
251   END SUBROUTINE lim_update2
252#else
253   !!----------------------------------------------------------------------
254   !!   Default option         Empty Module               No sea-ice model
255   !!----------------------------------------------------------------------
256CONTAINS
257   SUBROUTINE lim_update2     ! Empty routine
258   END SUBROUTINE lim_update2
259
260#endif
261
262END MODULE limupdate2
Note: See TracBrowser for help on using the repository browser.